The R Project SVN R-packages

Rev

Rev 1820 | Blame | Compare with Previous | Last modification | View Log | Download | RSS feed 

/* ========================================================================== */
/* === UMF_store_lu ========================================================= */
/* ========================================================================== */

/* -------------------------------------------------------------------------- */
/* 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                       */
/* -------------------------------------------------------------------------- */

/*
    Store the LU factors.  Called by the kernel.
    Returns TRUE if successful, FALSE if out of memory.
*/

#include "umf_internal.h"
#include "umf_mem_alloc_head_block.h"
#include "umf_get_memory.h"

/* ========================================================================== */

#ifdef DROP
GLOBAL Int UMF_store_lu_drop
#else
GLOBAL Int UMF_store_lu
#endif
(
    NumericType *Numeric,
    WorkType *Work
)
{
    /* ---------------------------------------------------------------------- */
    /* local variables */
    /* ---------------------------------------------------------------------- */

    Entry pivot_value ;
#ifdef DROP
    double droptol ;
#endif
    Entry *D, *Lval, *Uval, *Fl1, *Fl2, *Fu1, *Fu2,
    *Flublock, *Flblock, *Fublock ;
    Int i, k, fnr_curr, fnrows, fncols, row, col, pivrow, pivcol, *Frows,
    *Fcols, *Lpattern, *Upattern, *Lpos, *Upos, llen, ulen, fnc_curr, fnpiv,
    uilen, lnz, unz, nb, *Lilen,
    *Uilen, *Lip, *Uip, *Li, *Ui, pivcol_position, newLchain, newUchain,
    pivrow_position, p, size, lip, uip, lnzi, lnzx, unzx, lnz2i, lnz2x,
    unz2i, unz2x, zero_pivot, *Pivrow, *Pivcol, kk,
    Lnz [MAXNB] ;

#ifndef NDEBUG
    Int *Col_degree, *Row_degree ;
#endif

#ifdef DROP
    Int all_lnz, all_unz ;
    droptol = Numeric->droptol ;
#endif

    /* ---------------------------------------------------------------------- */
    /* get parameters */
    /* ---------------------------------------------------------------------- */

    fnrows = Work->fnrows ;
    fncols = Work->fncols ;
    fnpiv = Work->fnpiv ;

    Lpos = Numeric->Lpos ;
    Upos = Numeric->Upos ;
    Lilen = Numeric->Lilen ;
    Uilen = Numeric->Uilen ;

    Lip = Numeric->Lip ;
    Uip = Numeric->Uip ;
    D = Numeric->D ;

    Flublock = Work->Flublock ;
    Flblock  = Work->Flblock ;
    Fublock  = Work->Fublock ;

    fnr_curr = Work->fnr_curr ;
    fnc_curr = Work->fnc_curr ;
    Frows = Work->Frows ;
    Fcols = Work->Fcols ;

#ifndef NDEBUG
    Col_degree = Numeric->Cperm ;   /* for NON_PIVOTAL_COL macro */
    Row_degree = Numeric->Rperm ;   /* for NON_PIVOTAL_ROW macro */
#endif

    Lpattern = Work->Lpattern ;
    llen = Work->llen ;
    Upattern = Work->Upattern ;
    ulen = Work->ulen ;

    nb = Work->nb ;

#ifndef NDEBUG
    DEBUG1 (("\nSTORE LU: fnrows "ID
    " fncols "ID"\n", fnrows, fncols)) ;

    DEBUG2 (("\nFrontal matrix, including all space:\n"
        "fnr_curr "ID" fnc_curr "ID" nb    "ID"\n"
        "fnrows   "ID" fncols   "ID" fnpiv "ID"\n",
        fnr_curr, fnc_curr, nb, fnrows, fncols, fnpiv)) ;

    DEBUG2 (("\nJust the active part:\n")) ;
    DEBUG7 (("C  block: ")) ;
    UMF_dump_dense (Work->Fcblock,  fnr_curr, fnrows, fncols) ;
    DEBUG7 (("L  block: ")) ;
    UMF_dump_dense (Work->Flblock,  fnr_curr, fnrows, fnpiv);
    DEBUG7 (("U' block: ")) ;
    UMF_dump_dense (Work->Fublock,  fnc_curr, fncols, fnpiv) ;
    DEBUG7 (("LU block: ")) ;
    UMF_dump_dense (Work->Flublock, nb, fnpiv, fnpiv) ;
    DEBUG7 (("Current frontal matrix: (prior to store LU)\n")) ;
    UMF_dump_current_front (Numeric, Work, TRUE) ;
#endif

    Pivrow = Work->Pivrow ;
    Pivcol = Work->Pivcol ;

    /* ---------------------------------------------------------------------- */
    /* store the columns of L */
    /* ---------------------------------------------------------------------- */

    for (kk = 0 ; kk < fnpiv ; kk++)
    {

    /* ------------------------------------------------------------------ */
    /* one more pivot row and column is being stored into L and U */
    /* ------------------------------------------------------------------ */

    k = Work->npiv + kk ;

    /* ------------------------------------------------------------------ */
    /* find the kth pivot row and pivot column */
    /* ------------------------------------------------------------------ */

    pivrow = Pivrow [kk] ;
    pivcol = Pivcol [kk] ;

#ifndef NDEBUG
    ASSERT (pivrow >= 0 && pivrow < Work->n_row) ;
    ASSERT (pivcol >= 0 && pivcol < Work->n_col) ;

    DEBUGm4 ((
    "\n -------------------------------------------------------------"
    "Store LU: step " ID"\n", k))  ;
    ASSERT (k < MIN (Work->n_row, Work->n_col)) ;
    DEBUG2 (("Store column of L, k = "ID", llen "ID"\n", k, llen)) ;
    for (i = 0 ; i < llen ; i++)
    {
        row = Lpattern [i] ;
        ASSERT (row >= 0 && row < Work->n_row) ;
        DEBUG2 (("    Lpattern["ID"] "ID" Lpos "ID, i, row, Lpos [row])) ;
        if (row == pivrow) DEBUG2 ((" <- pivot row")) ;
        DEBUG2 (("\n")) ;
        ASSERT (i == Lpos [row]) ;
    }
#endif

    /* ------------------------------------------------------------------ */
    /* remove pivot row from L */
    /* ------------------------------------------------------------------ */

    /* remove pivot row index from current column of L */
    /* if a new Lchain starts, then all entries are removed later */
    DEBUG2 (("Removing pivrow from Lpattern, k = "ID"\n", k)) ;
    ASSERT (!NON_PIVOTAL_ROW (pivrow)) ;
    pivrow_position = Lpos [pivrow] ;
    if (pivrow_position != EMPTY)
    {
        /* place the last entry in the column in the */
        /* position of the pivot row index */
        ASSERT (pivrow == Lpattern [pivrow_position]) ;
        row = Lpattern [--llen] ;
        /* ASSERT (NON_PIVOTAL_ROW (row)) ; */
        Lpattern [pivrow_position] = row ;
        Lpos [row] = pivrow_position ;
        Lpos [pivrow] = EMPTY ;
    }

    /* ------------------------------------------------------------------ */
    /* store the pivot value, for the diagonal matrix D */
    /* ------------------------------------------------------------------ */

    /* kk-th column of LU block */
    Fl1 = Flublock + kk * nb ;

    /* kk-th column of L in the L block */
    Fl2 = Flblock + kk * fnr_curr ;

    /* kk-th pivot in frontal matrix located in Flublock [kk, kk] */
    pivot_value = Fl1 [kk] ;

    D [k] = pivot_value ;
    zero_pivot = IS_ZERO (pivot_value) ;

    DEBUG4 (("Pivot D["ID"]=", k)) ;
    EDEBUG4 (pivot_value) ;
    DEBUG4 (("\n")) ;

    /* ------------------------------------------------------------------ */
    /* count nonzeros in kth column of L */
    /* ------------------------------------------------------------------ */

    lnz = 0 ;
    lnz2i = 0 ;
    lnz2x = llen ;

#ifdef DROP
        all_lnz = 0 ;

        for (i = kk + 1 ; i < fnpiv ; i++)
        {
        Entry x ;
        double s ;
        x = Fl1 [i] ;
        if (IS_ZERO (x)) continue ;
        all_lnz++ ;
        APPROX_ABS (s, x) ;
        if (s <= droptol) continue ;
        lnz++ ;
        if (Lpos [Pivrow [i]] == EMPTY) lnz2i++ ;
        }

        for (i = 0 ; i < fnrows ; i++)
        {
        Entry x ;
        double s ;
        x = Fl2 [i] ;
        if (IS_ZERO (x)) continue ;
        all_lnz++ ;
        APPROX_ABS (s, x) ;
        if (s <= droptol) continue ;
        lnz++ ;
        if (Lpos [Frows [i]] == EMPTY) lnz2i++ ;
        }

#else

        for (i = kk + 1 ; i < fnpiv ; i++)
        {
        if (IS_ZERO (Fl1 [i])) continue ;
        lnz++ ;
        if (Lpos [Pivrow [i]] == EMPTY) lnz2i++ ;
        }

        for (i = 0 ; i < fnrows ; i++)
        {
        if (IS_ZERO (Fl2 [i])) continue ;
        lnz++ ;
        if (Lpos [Frows [i]] == EMPTY) lnz2i++ ;
        }

#endif

    lnz2x += lnz2i ;

    /* determine if we start a new Lchain or continue the old one */
    if (llen == 0 || zero_pivot)
    {
        /* llen == 0 means there is no prior Lchain */
        /* D [k] == 0 means the pivot column is empty */
        newLchain = TRUE ;
    }
    else
    {
        newLchain =
            /* storage for starting a new Lchain */
            UNITS (Entry, lnz) + UNITS (Int, lnz)
        <=
            /* storage for continuing a prior Lchain */
            UNITS (Entry, lnz2x) + UNITS (Int, lnz2i) ;
    }

    if (newLchain)
    {
        /* start a new chain for column k of L */
        DEBUG2 (("Start new Lchain, k = "ID"\n", k)) ;

        pivrow_position = EMPTY ;

        /* clear the prior Lpattern */
        for (i = 0 ; i < llen ; i++)
        {
        row = Lpattern [i] ;
        Lpos [row] = EMPTY ;
        }
        llen = 0 ;

        lnzi = lnz ;
        lnzx = lnz ;
    }
    else
    {
        /* continue the prior Lchain */
        DEBUG2 (("Continue  Lchain, k = "ID"\n", k)) ;
        lnzi = lnz2i ;
        lnzx = lnz2x ;
    }

    /* ------------------------------------------------------------------ */
    /* allocate space for the column of L */
    /* ------------------------------------------------------------------ */

    size = UNITS (Int, lnzi) + UNITS (Entry, lnzx) ;

#ifndef NDEBUG
    UMF_allocfail = FALSE ;
    if (UMF_gprob > 0)
    {
        double rrr = ((double) (rand ( ))) / (((double) RAND_MAX) + 1) ;
        DEBUG4 (("Check random %e %e\n", rrr, UMF_gprob)) ;
        UMF_allocfail = rrr < UMF_gprob ;
        if (UMF_allocfail) DEBUGm2 (("Random garbage coll. (store LU)\n"));
    }
#endif

    p = UMF_mem_alloc_head_block (Numeric, size) ;
    if (!p)
    {
        Int r2, c2 ;
        /* Do garbage collection, realloc, and try again. */
        /* Note that there are pivot rows/columns in current front. */
        if (Work->do_grow)
        {
        /* full compaction of current frontal matrix, since
         * UMF_grow_front will be called next anyway. */
        r2 = fnrows ;
        c2 = fncols ;
        }
        else
        {
        /* partial compaction. */
        r2 = MAX (fnrows, Work->fnrows_new + 1) ;
        c2 = MAX (fncols, Work->fncols_new + 1) ;
        }
        DEBUGm3 (("get_memory from umf_store_lu:\n")) ;
        if (!UMF_get_memory (Numeric, Work, size, r2, c2, TRUE))
        {
        DEBUGm4 (("out of memory: store LU (1)\n")) ;
        return (FALSE) ;    /* out of memory */
        }
        p = UMF_mem_alloc_head_block (Numeric, size) ;
        if (!p)
        {
        DEBUGm4 (("out of memory: store LU (2)\n")) ;
        return (FALSE) ;    /* out of memory */
        }
        /* garbage collection may have moved the current front */
        fnc_curr = Work->fnc_curr ;
        fnr_curr = Work->fnr_curr ;
        Flublock = Work->Flublock ;
        Flblock  = Work->Flblock ;
        Fublock  = Work->Fublock ;
        Fl1 = Flublock + kk * nb ;
        Fl2 = Flblock  + kk * fnr_curr ;
    }

    /* ------------------------------------------------------------------ */
    /* store the column of L */
    /* ------------------------------------------------------------------ */

    lip = p ;

    Li = (Int *) (Numeric->Memory + p) ;
    p += UNITS (Int, lnzi) ;
    Lval = (Entry *) (Numeric->Memory + p) ;
    p += UNITS (Entry, lnzx) ;

    for (i = 0 ; i < lnzx ; i++)
    {
        CLEAR (Lval [i]) ;
    }

    /* store the numerical entries */

    if (newLchain)
    {
        /* flag the first column in the Lchain by negating Lip [k] */
        lip = -lip ;

        ASSERT (llen == 0) ;

#ifdef DROP

        for (i = kk + 1 ; i < fnpiv ; i++)
        {
            Entry x ;
            double s ;
            Int row2, pos ;
            x = Fl1 [i] ;
            APPROX_ABS (s, x) ;
            if (s <= droptol) continue ;
            row2 = Pivrow [i] ;
            pos = llen++ ;
            Lpattern [pos] = row2 ;
            Lpos [row2] = pos ;
            Li [pos] = row2 ;
            Lval [pos] = x ;
        }

        for (i = 0 ; i < fnrows ; i++)
        {
            Entry x ;
            double s ;
            Int row2, pos ;
            x = Fl2 [i] ;
            APPROX_ABS (s, x) ;
            if (s <= droptol) continue ;
            row2 = Frows [i] ;
            pos = llen++ ;
            Lpattern [pos] = row2 ;
            Lpos [row2] = pos ;
            Li [pos] = row2 ;
            Lval [pos] = x ;
        }

#else

        for (i = kk + 1 ; i < fnpiv ; i++)
        {
            Entry x ;
            Int row2, pos ;
            x = Fl1 [i] ;
            if (IS_ZERO (x)) continue ;
            row2 = Pivrow [i] ;
            pos = llen++ ;
            Lpattern [pos] = row2 ;
            Lpos [row2] = pos ;
            Li [pos] = row2 ;
            Lval [pos] = x ;
        }

        for (i = 0 ; i < fnrows ; i++)
        {
            Entry x ;
            Int row2, pos ;
            x = Fl2 [i] ;
            if (IS_ZERO (x)) continue ;
            row2 = Frows [i] ;
            pos = llen++ ;
            Lpattern [pos] = row2 ;
            Lpos [row2] = pos ;
            Li [pos] = row2 ;
            Lval [pos] = x ;
        }

#endif

    }
    else
    {
        ASSERT (llen > 0) ;

#ifdef DROP

        for (i = kk + 1 ; i < fnpiv ; i++)
        {
            Entry x ;
            double s ;
            Int row2, pos ;
            x = Fl1 [i] ;
            APPROX_ABS (s, x) ;
            if (s <= droptol) continue ;
            row2 = Pivrow [i] ;
            pos = Lpos [row2] ;
            if (pos == EMPTY)
            {
            pos = llen++ ;
            Lpattern [pos] = row2 ;
            Lpos [row2] = pos ;
            *Li++ = row2 ;
            }
            Lval [pos] = x ;
        }

        for (i = 0 ; i < fnrows ; i++)
        {
            Entry x ;
            double s ;
            Int row2, pos ;
            x = Fl2 [i] ;
            APPROX_ABS (s, x) ;
            if (s <= droptol) continue ;
            row2 = Frows [i] ;
            pos = Lpos [row2] ;
            if (pos == EMPTY)
            {
            pos = llen++ ;
            Lpattern [pos] = row2 ;
            Lpos [row2] = pos ;
            *Li++ = row2 ;
            }
            Lval [pos] = x ;
        }

#else

        for (i = kk + 1 ; i < fnpiv ; i++)
        {
            Entry x ;
            Int row2, pos ;
            x = Fl1 [i] ;
            if (IS_ZERO (x)) continue ;
            row2 = Pivrow [i] ;
            pos = Lpos [row2] ;
            if (pos == EMPTY)
            {
            pos = llen++ ;
            Lpattern [pos] = row2 ;
            Lpos [row2] = pos ;
            *Li++ = row2 ;
            }
            Lval [pos] = x ;
        }

        for (i = 0 ; i < fnrows ; i++)
        {
            Entry x ;
            Int row2, pos ;
            x = Fl2 [i] ;
            if (IS_ZERO (x)) continue ;
            row2 = Frows [i] ;
            pos = Lpos [row2] ;
            if (pos == EMPTY)
            {
            pos = llen++ ;
            Lpattern [pos] = row2 ;
            Lpos [row2] = pos ;
            *Li++ = row2 ;
            }
            Lval [pos] = x ;
        }

#endif

    }
    DEBUG4 (("llen "ID" lnzx "ID"\n", llen, lnzx)) ;
    ASSERT (llen == lnzx) ;
    ASSERT (lnz <= llen) ;
    DEBUG4 (("lnz "ID" \n", lnz)) ;

#ifdef DROP

        DEBUG4 (("all_lnz "ID" \n", all_lnz)) ;
        ASSERT (lnz <= all_lnz) ;
        Numeric->lnz += lnz ;
        Numeric->all_lnz += all_lnz ;
        Lnz [kk] = all_lnz ;

#else

        Numeric->lnz += lnz ;
        Numeric->all_lnz += lnz ;
        Lnz [kk] = lnz ;
#endif

    Numeric->nLentries += lnzx ;
    Work->llen = llen ;
    Numeric->isize += lnzi ;

    /* ------------------------------------------------------------------ */
    /* the pivot column is fully assembled and scaled, and is now the */
    /* k-th column of L */
    /* ------------------------------------------------------------------ */

    Lpos [pivrow] = pivrow_position ;   /* not aliased */
    Lip [pivcol] = lip ;            /* aliased with Col_tuples */
    Lilen [pivcol] = lnzi ;         /* aliased with Col_tlen */

    }

    /* ---------------------------------------------------------------------- */
    /* store the rows of U */
    /* ---------------------------------------------------------------------- */

    for (kk = 0 ; kk < fnpiv ; kk++)
    {

    /* ------------------------------------------------------------------ */
    /* one more pivot row and column is being stored into L and U */
    /* ------------------------------------------------------------------ */

    k = Work->npiv + kk ;

    /* ------------------------------------------------------------------ */
    /* find the kth pivot row and pivot column */
    /* ------------------------------------------------------------------ */

    pivrow = Pivrow [kk] ;
    pivcol = Pivcol [kk] ;

#ifndef NDEBUG
    ASSERT (pivrow >= 0 && pivrow < Work->n_row) ;
    ASSERT (pivcol >= 0 && pivcol < Work->n_col) ;

    DEBUG2 (("Store row of U, k = "ID", ulen "ID"\n", k, ulen)) ;
    for (i = 0 ; i < ulen ; i++)
    {
        col = Upattern [i] ;
        DEBUG2 (("    Upattern["ID"] "ID, i, col)) ;
        if (col == pivcol) DEBUG2 ((" <- pivot col")) ;
        DEBUG2 (("\n")) ;
        ASSERT (col >= 0 && col < Work->n_col) ;
        ASSERT (i == Upos [col]) ;
    }
#endif

    /* ------------------------------------------------------------------ */
    /* get the pivot value, for the diagonal matrix D */
    /* ------------------------------------------------------------------ */

    zero_pivot = IS_ZERO (D [k]) ;

    /* ------------------------------------------------------------------ */
    /* count the nonzeros in the row of U */
    /* ------------------------------------------------------------------ */

    /* kk-th row of U in the LU block */
    Fu1 = Flublock + kk ;

    /* kk-th row of U in the U block */
    Fu2 = Fublock + kk * fnc_curr ;

    unz = 0 ;
    unz2i = 0 ;
    unz2x = ulen ;
    DEBUG2 (("unz2x is "ID", lnzx "ID"\n", unz2x, lnzx)) ;

    /* if row k does not end a Uchain, pivcol not included in ulen */
    ASSERT (!NON_PIVOTAL_COL (pivcol)) ;
    pivcol_position = Upos [pivcol] ;
    if (pivcol_position != EMPTY)
    {
        unz2x-- ;
        DEBUG2 (("(exclude pivcol) unz2x is now "ID"\n", unz2x)) ;
    }

    ASSERT (unz2x >= 0) ;

#ifdef DROP
        all_unz = 0 ;

        for (i = kk + 1 ; i < fnpiv ; i++)
        {
        Entry x ;
        double s ;
        x = Fu1 [i*nb] ;
        if (IS_ZERO (x)) continue ;
        all_unz++ ;
        APPROX_ABS (s, x) ;
        if (s <= droptol) continue ;
        unz++ ;
        if (Upos [Pivcol [i]] == EMPTY) unz2i++ ;
        }

        for (i = 0 ; i < fncols ; i++)
        {
        Entry x ;
        double s ;
        x = Fu2 [i] ;
        if (IS_ZERO (x)) continue ;
        all_unz++ ;
        APPROX_ABS (s, x) ;
        if (s <= droptol) continue ;
        unz++ ;
        if (Upos [Fcols [i]] == EMPTY) unz2i++ ;
        }

#else

        for (i = kk + 1 ; i < fnpiv ; i++)
        {
        if (IS_ZERO (Fu1 [i*nb])) continue ;
        unz++ ;
        if (Upos [Pivcol [i]] == EMPTY) unz2i++ ;
        }

        for (i = 0 ; i < fncols ; i++)
        {
        if (IS_ZERO (Fu2 [i])) continue ;
        unz++ ;
        if (Upos [Fcols [i]] == EMPTY) unz2i++ ;
        }

#endif

    unz2x += unz2i ;

    ASSERT (IMPLIES (k == 0, ulen == 0)) ;

    /* determine if we start a new Uchain or continue the old one */
    if (ulen == 0 || zero_pivot)
    {
        /* ulen == 0 means there is no prior Uchain */
        /* D [k] == 0 means the matrix is singular (pivot row might */
        /* not be empty, however, but start a new Uchain to prune zero */
        /* entries for the deg > 0 test in UMF_u*solve) */
        newUchain = TRUE ;
    }
    else
    {
        newUchain =
            /* approximate storage for starting a new Uchain */
            UNITS (Entry, unz) + UNITS (Int, unz)
        <=
            /* approximate storage for continuing a prior Uchain */
            UNITS (Entry, unz2x) + UNITS (Int, unz2i) ;

        /* this would be exact, except for the Int to Unit rounding, */
        /* because the Upattern is stored only at the end of the Uchain */
    }

    /* ------------------------------------------------------------------ */
    /* allocate space for the row of U */
    /* ------------------------------------------------------------------ */

    size = 0 ;
    if (newUchain)
    {
        /* store the pattern of the last row in the prior Uchain */
        size += UNITS (Int, ulen) ;
        unzx = unz ;
    }
    else
    {
        unzx = unz2x ;
    }
    size += UNITS (Entry, unzx) ;

#ifndef NDEBUG
    UMF_allocfail = FALSE ;
    if (UMF_gprob > 0)
    {
        double rrr = ((double) (rand ( ))) / (((double) RAND_MAX) + 1) ;
        DEBUG4 (("Check random %e %e\n", rrr, UMF_gprob)) ;
        UMF_allocfail = rrr < UMF_gprob ;
        if (UMF_allocfail) DEBUGm2 (("Random garbage coll. (store LU)\n"));
    }
#endif

    p = UMF_mem_alloc_head_block (Numeric, size) ;
    if (!p)
    {
        Int r2, c2 ;
        /* Do garbage collection, realloc, and try again. */
        /* Note that there are pivot rows/columns in current front. */
        if (Work->do_grow)
        {
        /* full compaction of current frontal matrix, since
         * UMF_grow_front will be called next anyway. */
        r2 = fnrows ;
        c2 = fncols ;
        }
        else
        {
        /* partial compaction. */
        r2 = MAX (fnrows, Work->fnrows_new + 1) ;
        c2 = MAX (fncols, Work->fncols_new + 1) ;
        }
        DEBUGm3 (("get_memory from umf_store_lu:\n")) ;
        if (!UMF_get_memory (Numeric, Work, size, r2, c2, TRUE))
        {
        /* :: get memory, column of L :: */
        DEBUGm4 (("out of memory: store LU (1)\n")) ;
        return (FALSE) ;    /* out of memory */
        }
        p = UMF_mem_alloc_head_block (Numeric, size) ;
        if (!p)
        {
        /* :: out of memory, column of U :: */
        DEBUGm4 (("out of memory: store LU (2)\n")) ;
        return (FALSE) ;    /* out of memory */
        }
        /* garbage collection may have moved the current front */
        fnc_curr = Work->fnc_curr ;
        fnr_curr = Work->fnr_curr ;
        Flublock = Work->Flublock ;
        Flblock  = Work->Flblock ;
        Fublock  = Work->Fublock ;
        Fu1 = Flublock + kk ;
        Fu2 = Fublock  + kk * fnc_curr ;
    }

    /* ------------------------------------------------------------------ */
    /* store the row of U */
    /* ------------------------------------------------------------------ */

    uip = p ;

    if (newUchain)
    {
        /* starting a new Uchain - flag this by negating Uip [k] */
        uip = -uip ;
        DEBUG2 (("Start new Uchain, k = "ID"\n", k)) ;

        pivcol_position = EMPTY ;

        /* end the prior Uchain */
        /* save the current Upattern, and then */
        /* clear it and start a new Upattern */
        DEBUG2 (("Ending prior chain, k-1 = "ID"\n", k-1)) ;
        uilen = ulen ;
        Ui = (Int *) (Numeric->Memory + p) ;
        Numeric->isize += ulen ;
        p += UNITS (Int, ulen) ;
        for (i = 0 ; i < ulen ; i++)
        {
        col = Upattern [i] ;
        ASSERT (col >= 0 && col < Work->n_col) ;
        Upos [col] = EMPTY ;
        Ui [i] = col ;
        }

        ulen = 0 ;

    }
    else
    {
        /* continue the prior Uchain */
        DEBUG2 (("Continue  Uchain, k = "ID"\n", k)) ;
        ASSERT (k > 0) ;

        /* remove pivot col index from current row of U */
        /* if a new Uchain starts, then all entries are removed later */
        DEBUG2 (("Removing pivcol from Upattern, k = "ID"\n", k)) ;

        if (pivcol_position != EMPTY)
        {
        /* place the last entry in the row in the */
        /* position of the pivot col index */
        ASSERT (pivcol == Upattern [pivcol_position]) ;
        col = Upattern [--ulen] ;
        ASSERT (col >= 0 && col < Work->n_col) ;
        Upattern [pivcol_position] = col ;
        Upos [col] = pivcol_position ;
        Upos [pivcol] = EMPTY ;
        }

        /* this row continues the Uchain.  Keep track of how much */
        /* to trim from the k-th length to get the length of the */
        /* (k-1)st row of U */
        uilen = unz2i ;

    }

    Uval = (Entry *) (Numeric->Memory + p) ;
    /* p += UNITS (Entry, unzx), no need to increment p */

    for (i = 0 ; i < unzx ; i++)
    {
        CLEAR (Uval [i]) ;
    }

    if (newUchain)
    {
        ASSERT (ulen == 0) ;

#ifdef DROP

        for (i = kk + 1 ; i < fnpiv ; i++)
        {
            Entry x ;
            double s ;
            Int col2, pos ;
            x = Fu1 [i*nb] ;
            APPROX_ABS (s, x) ;
            if (s <= droptol) continue ;
            col2 = Pivcol [i] ;
            pos = ulen++ ;
            Upattern [pos] = col2 ;
            Upos [col2] = pos ;
            Uval [pos] = x ;
        }

        for (i = 0 ; i < fncols ; i++)
        {
            Entry x ;
            double s ;
            Int col2, pos ;
            x = Fu2 [i] ;
            APPROX_ABS (s, x) ;
            if (s <= droptol) continue ;
            col2 = Fcols [i] ;
            pos = ulen++ ;
            Upattern [pos] = col2 ;
            Upos [col2] = pos ;
            Uval [pos] = x ;
        }

#else

        for (i = kk + 1 ; i < fnpiv ; i++)
        {
            Entry x ;
            Int col2, pos ;
            x = Fu1 [i*nb] ;
            if (IS_ZERO (x)) continue ;
            col2 = Pivcol [i] ;
            pos = ulen++ ;
            Upattern [pos] = col2 ;
            Upos [col2] = pos ;
            Uval [pos] = x ;
        }

        for (i = 0 ; i < fncols ; i++)
        {
            Entry x ;
            Int col2, pos ;
            x = Fu2 [i] ;
            if (IS_ZERO (x)) continue ;
            col2 = Fcols [i] ;
            pos = ulen++ ;
            Upattern [pos] = col2 ;
            Upos [col2] = pos ;
            Uval [pos] = x ;
        }

#endif

    }
    else
    {

        ASSERT (ulen > 0) ;

        /* store the numerical entries and find new nonzeros */

#ifdef DROP

        for (i = kk + 1 ; i < fnpiv ; i++)
        {
            Entry x ;
            double s ;
            Int col2, pos ;
            x = Fu1 [i*nb] ;
            APPROX_ABS (s, x) ;
            if (s <= droptol) continue ;
            col2 = Pivcol [i] ;
            pos = Upos [col2] ;
            if (pos == EMPTY)
            {
            pos = ulen++ ;
            Upattern [pos] = col2 ;
            Upos [col2] = pos ;
            }
            Uval [pos] = x ;
        }

        for (i = 0 ; i < fncols ; i++)
        {
            Entry x ;
            double s ;
            Int col2, pos ;
            x = Fu2 [i] ;
            APPROX_ABS (s, x) ;
            if (s <= droptol) continue ;
            col2 = Fcols [i] ;
            pos = Upos [col2] ;
            if (pos == EMPTY)
            {
            pos = ulen++ ;
            Upattern [pos] = col2 ;
            Upos [col2] = pos ;
            }
            Uval [pos] = x ;
        }

#else

        for (i = kk + 1 ; i < fnpiv ; i++)
        {
            Entry x ;
            Int col2, pos ;
            x = Fu1 [i*nb] ;
            if (IS_ZERO (x)) continue ;
            col2 = Pivcol [i] ;
            pos = Upos [col2] ;
            if (pos == EMPTY)
            {
            pos = ulen++ ;
            Upattern [pos] = col2 ;
            Upos [col2] = pos ;
            }
            Uval [pos] = x ;
        }

        for (i = 0 ; i < fncols ; i++)
        {
            Entry x ;
            Int col2, pos ;
            x = Fu2 [i] ;
            if (IS_ZERO (x)) continue ;
            col2 = Fcols [i] ;
            pos = Upos [col2] ;
            if (pos == EMPTY)
            {
            pos = ulen++ ;
            Upattern [pos] = col2 ;
            Upos [col2] = pos ;
            }
            Uval [pos] = x ;
        }

#endif

    }

    ASSERT (ulen == unzx) ;
    ASSERT (unz <= ulen) ;
    DEBUG4 (("unz "ID" \n", unz)) ;

#ifdef DROP

        DEBUG4 (("all_unz "ID" \n", all_unz)) ;
        ASSERT (unz <= all_unz) ;
        Numeric->unz += unz ;
        Numeric->all_unz += all_unz ;
        /* count the "true" flops, based on LU pattern only */
        Numeric->flops += DIV_FLOPS * Lnz [kk]  /* scale pivot column */
        + MULTSUB_FLOPS * (Lnz [kk] * all_unz) ;    /* outer product */

#else

        Numeric->unz += unz ;
        Numeric->all_unz += unz ;
        /* count the "true" flops, based on LU pattern only */
        Numeric->flops += DIV_FLOPS * Lnz [kk]  /* scale pivot column */
        + MULTSUB_FLOPS * (Lnz [kk] * unz) ;    /* outer product */
#endif

    Numeric->nUentries += unzx ;
    Work->ulen = ulen ;
    DEBUG1 (("Work->ulen = "ID" at end of pivot step, k: "ID"\n", ulen, k));

    /* ------------------------------------------------------------------ */
    /* the pivot row is the k-th row of U */
    /* ------------------------------------------------------------------ */

    Upos [pivcol] = pivcol_position ;   /* not aliased */
    Uip [pivrow] = uip ;            /* aliased with Row_tuples */
    Uilen [pivrow] = uilen ;        /* aliased with Row_tlen */

    }

    /* ---------------------------------------------------------------------- */
    /* no more pivots in frontal working array */
    /* ---------------------------------------------------------------------- */

    Work->npiv += fnpiv ;
    Work->fnpiv = 0 ;
    Work->fnzeros = 0 ;
    return (TRUE) ;
}