The R Project SVN R-packages

/* ========================================================================== */
/* === UMF_assemble ========================================================= */
/* ========================================================================== */

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

/*  Degree update and numerical assembly.  This is compiled twice (with and
 *  without FIXQ) for each real/complex int/long version, for a total of 8
 *  versions.*/

#include "umf_internal.h"
#include "umf_mem_free_tail_block.h"

/* ========================================================================== */
/* === row_assemble ========================================================= */
/* ========================================================================== */

PRIVATE void row_assemble
(
    Int row,
    NumericType *Numeric,
    WorkType *Work
)
{

    Entry *S, *Fcblock, *Frow ;
    Int tpi, e, *E, *Fcpos, *Frpos, *Row_degree, *Row_tuples, *Row_tlen, rdeg0,
    f, nrows, ncols, *Rows, *Cols, col, ncolsleft, j ;
    Tuple *tp, *tp1, *tp2, *tpend ;
    Unit *Memory, *p ;
    Element *ep ;

#ifndef FIXQ
    Int *Col_degree ;
    Col_degree = Numeric->Cperm ;
#endif

    Row_tuples = Numeric->Uip ;
    tpi = Row_tuples [row] ;
    if (!tpi) return ;

    Memory = Numeric->Memory ;
    E = Work->E ;
    Fcpos = Work->Fcpos ;
    Frpos = Work->Frpos ;
    Row_degree = Numeric->Rperm ;
    Row_tlen   = Numeric->Uilen ;
    E = Work->E ;
    Memory = Numeric->Memory ;
    rdeg0 = Work->rdeg0 ;
    Fcblock = Work->Fcblock ;

#ifndef NDEBUG
    DEBUG6 (("SCAN2-row: "ID"\n", row)) ;
    UMF_dump_rowcol (0, Numeric, Work, row, FALSE) ;
#endif

    ASSERT (NON_PIVOTAL_ROW (row)) ;

    tp = (Tuple *) (Memory + tpi) ;
    tp1 = tp ;
    tp2 = tp ;
    tpend = tp + Row_tlen [row] ;
    for ( ; tp < tpend ; tp++)
    {
    e = tp->e ;
    ASSERT (e > 0 && e <= Work->nel) ;
    if (!E [e]) continue ;  /* element already deallocated */
    f = tp->f ;
    p = Memory + E [e] ;
    ep = (Element *) p ;
    p += UNITS (Element, 1) ;
    Cols = (Int *) p ;
    Rows = Cols + ep->ncols ;
    if (Rows [f] == EMPTY) continue ;   /* row already assembled */
    ASSERT (row == Rows [f] && row >= 0 && row < Work->n_row) ;

    if (ep->rdeg == rdeg0)
    {
        /* ------------------------------------------------------ */
        /* this is an old Lson - assemble just one row */
        /* ------------------------------------------------------ */

        /* flag the row as assembled from the Lson */
        Rows [f] = EMPTY ;

        nrows = ep->nrows ;
        ncols = ep->ncols ;

        p += UNITS (Int, ncols + nrows) ;
        S = ((Entry *) p) + f ;

        DEBUG6 (("Old LSON: "ID"\n", e)) ;
#ifndef NDEBUG
        UMF_dump_element (Numeric, Work, e, FALSE) ;
#endif

        ncolsleft = ep->ncolsleft ;

        Frow = Fcblock + Frpos [row] ;
        DEBUG6 (("LSON found (in scan2-row): "ID"\n", e)) ;

        Row_degree [row] -= ncolsleft ;

        if (ncols == ncolsleft)
        {
        /* -------------------------------------------------- */
        /* no columns assembled out this Lson yet */
        /* -------------------------------------------------- */

#pragma ivdep
        for (j = 0 ; j < ncols ; j++)
        {
            col = Cols [j] ;
            ASSERT (col >= 0 && col < Work->n_col) ;
#ifndef FIXQ
            Col_degree [col] -- ;
#endif
            /* Frow [Fcpos [col]] += *S ; */
            ASSEMBLE (Frow [Fcpos [col]], *S) ;
            S += nrows ;
        }

        }
        else
        {
        /* -------------------------------------------------- */
        /* some columns have been assembled out of this Lson */
        /* -------------------------------------------------- */

#pragma ivdep
        for (j = 0 ; j < ncols ; j++)
        {
            col = Cols [j] ;
            if (col >= 0)
            {
            ASSERT (col < Work->n_col) ;
#ifndef FIXQ
            Col_degree [col] -- ;
#endif
            /* Frow [Fcpos [col]] += *S ; */
            ASSEMBLE (Frow [Fcpos [col]], *S) ;
            }
            S += nrows ;
        }

        }
        ep->nrowsleft-- ;
        ASSERT (ep->nrowsleft > 0) ;
    }
    else
    {
        *tp2++ = *tp ;  /* leave the tuple in the list */
    }
    }
    Row_tlen [row] = tp2 - tp1 ;

#ifndef NDEBUG
    DEBUG7 (("row assembled in scan2-row: "ID"\n", row)) ;
    UMF_dump_rowcol (0, Numeric, Work, row, FALSE) ;
    DEBUG7 (("Current frontal matrix: (scan 1b)\n")) ;
    UMF_dump_current_front (Numeric, Work, TRUE) ;
#endif
}

/* ========================================================================== */
/* === col_assemble ========================================================= */
/* ========================================================================== */

PRIVATE void col_assemble
(
    Int col,
    NumericType *Numeric,
    WorkType *Work
)
{

    Entry *S, *Fcblock, *Fcol ;
    Int tpi, e, *E, *Fcpos, *Frpos, *Row_degree, *Col_tuples, *Col_tlen, cdeg0,
    f, nrows, ncols, *Rows, *Cols, row, nrowsleft, i ;
    Tuple *tp, *tp1, *tp2, *tpend ;
    Unit *Memory, *p ;
    Element *ep ;

#if !defined (FIXQ) || !defined (NDEBUG)
    Int *Col_degree ;
    Col_degree = Numeric->Cperm ;
#endif

    Col_tuples = Numeric->Lip ;
    tpi = Col_tuples [col] ;
    if (!tpi) return ;

    Memory = Numeric->Memory ;
    E = Work->E ;
    Fcpos = Work->Fcpos ;
    Frpos = Work->Frpos ;
    Row_degree = Numeric->Rperm ;
    Col_tlen   = Numeric->Lilen ;
    E = Work->E ;
    Memory = Numeric->Memory ;
    cdeg0 = Work->cdeg0 ;
    Fcblock = Work->Fcblock ;

    DEBUG6 (("SCAN2-col: "ID"\n", col)) ;
#ifndef NDEBUG
    UMF_dump_rowcol (1, Numeric, Work, col, FALSE) ;
#endif

    ASSERT (NON_PIVOTAL_COL (col)) ;
    tp = (Tuple *) (Memory + tpi) ;
    tp1 = tp ;
    tp2 = tp ;
    tpend = tp + Col_tlen [col] ;
    for ( ; tp < tpend ; tp++)
    {
    e = tp->e ;
    ASSERT (e > 0 && e <= Work->nel) ;
    if (!E [e]) continue ;  /* element already deallocated */
    f = tp->f ;
    p = Memory + E [e] ;
    ep = (Element *) p ;
    p += UNITS (Element, 1) ;
    Cols = (Int *) p ;

    if (Cols [f] == EMPTY) continue ;   /* col already assembled */
    ASSERT (col == Cols [f] && col >= 0 && col < Work->n_col) ;

    if (ep->cdeg == cdeg0)
    {
        /* ------------------------------------------------------ */
        /* this is an old Uson - assemble just one col */
        /* ------------------------------------------------------ */

        /* flag the col as assembled from the Uson */
        Cols [f] = EMPTY ;

        nrows = ep->nrows ;
        ncols = ep->ncols ;
        Rows = Cols + ncols ;
        p += UNITS (Int, ncols + nrows) ;
        S = ((Entry *) p) + f * nrows ;

        DEBUG6 (("Old USON: "ID"\n", e)) ;
#ifndef NDEBUG
        UMF_dump_element (Numeric, Work, e, FALSE) ;
#endif

        nrowsleft = ep->nrowsleft ;

        Fcol = Fcblock + Fcpos [col] ;
        DEBUG6 (("USON found (in scan2-col): "ID"\n", e)) ;
#ifndef FIXQ
        Col_degree [col] -= nrowsleft ;
#endif
        if (nrows == nrowsleft)
        {
        /* -------------------------------------------------- */
        /* no rows assembled out of this Uson yet */
        /* -------------------------------------------------- */

#pragma ivdep
        for (i = 0 ; i < nrows ; i++)
        {
            row = Rows [i] ;
            ASSERT (row >= 0 && row < Work->n_row) ;
            Row_degree [row]-- ;
            /* Fcol [Frpos [row]] += S [i] ; */
            ASSEMBLE (Fcol [Frpos [row]], S [i]) ;
        }
        }
        else
        {
        /* -------------------------------------------------- */
        /* some rows have been assembled out of this Uson */
        /* -------------------------------------------------- */

#pragma ivdep
        for (i = 0 ; i < nrows ; i++)
        {
            row = Rows [i] ;
            if (row >= 0)
            {
            ASSERT (row < Work->n_row) ;
            Row_degree [row]-- ;
            /* Fcol [Frpos [row]] += S [i] ; */
            ASSEMBLE (Fcol [Frpos [row]], S [i]) ;
            }
        }
        }
        ep->ncolsleft-- ;
        ASSERT (ep->ncolsleft > 0) ;
    }
    else
    {
        *tp2++ = *tp ;  /* leave the tuple in the list */
    }
    }
    Col_tlen [col] = tp2 - tp1 ;

#ifndef NDEBUG
    DEBUG7 (("Column assembled in scan2-col: "ID"\n", col)) ;
    UMF_dump_rowcol (1, Numeric, Work, col, FALSE) ;
    DEBUG7 (("Current frontal matrix: after scan2-col\n")) ;
    UMF_dump_current_front (Numeric, Work, TRUE) ;
#endif

}


/* ========================================================================== */
/* === UMF_assemble / UMF_assemble_fixq ===================================== */
/* ========================================================================== */

#ifndef FIXQ
GLOBAL void UMF_assemble
#else
GLOBAL void UMF_assemble_fixq
#endif
(
    NumericType *Numeric,
    WorkType *Work
)
{
    /* ---------------------------------------------------------------------- */
    /* local variables */
    /* ---------------------------------------------------------------------- */

    Int e, i, row, col, i2, nrows, ncols, f, tpi, extcdeg, extrdeg, rdeg0,
    cdeg0, son_list, next, nrows_to_assemble,
    ncols_to_assemble, ngetrows, j, j2,
    nrowsleft,  /* number of rows remaining in S */
    ncolsleft,  /* number of columns remaining in S */
    prior_Lson, prior_Uson, *E, *Cols, *Rows, *Wm, *Woo,
    *Row_tuples, *Row_degree, *Row_tlen,
    *Col_tuples, *Col_tlen ;
    Unit *Memory, *p ;
    Element *ep ;
    Tuple *tp, *tp1, *tp2, *tpend ;
    Entry
    *S,     /* a pointer into the contribution block of a son */
    *Fcblock,   /* current contribution block */
    *Fcol ;     /* a column of FC */
    Int *Frpos,
    *Fcpos,
    fnrows,     /* number of rows in contribution block in F */
    fncols ;    /* number of columns in contribution block in F */

#if !defined (FIXQ) || !defined (NDEBUG)
    Int *Col_degree ;
#endif

#ifndef NDEBUG
    Int n_row, n_col ;
    n_row = Work->n_row ;
    n_col = Work->n_col ;
    DEBUG3 (("::Assemble SCANS 1-4\n")) ;
    UMF_dump_current_front (Numeric, Work, TRUE) ;
#endif

#if !defined (FIXQ) || !defined (NDEBUG)
    Col_degree = Numeric->Cperm ;   /* not updated if FIXQ is true */
#endif

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

    fncols = Work->fncols ;
    fnrows = Work->fnrows ;
    Fcpos = Work->Fcpos ;
    Frpos = Work->Frpos ;
    Row_degree = Numeric->Rperm ;
    Row_tuples = Numeric->Uip ;
    Row_tlen   = Numeric->Uilen ;
    Col_tuples = Numeric->Lip ;
    Col_tlen   = Numeric->Lilen ;
    E = Work->E ;
    Memory = Numeric->Memory ;
    Wm = Work->Wm ;
    Woo = Work->Woo ;
    rdeg0 = Work->rdeg0 ;
    cdeg0 = Work->cdeg0 ;

#ifndef NDEBUG
    DEBUG6 (("============================================\n")) ;
    DEBUG6 (("Degree update, assembly.\n")) ;
    DEBUG6 (("pivot row pattern:  fncols="ID"\n", fncols)) ;
    for (j = 0 ; j < fncols ; j++)
    {
    col = Work->Fcols [j] ;
    DEBUG6 ((ID" ", col)) ;
    ASSERT (Fcpos [col] == j * Work->fnr_curr) ;
    ASSERT (NON_PIVOTAL_COL (col)) ;
    }
    ASSERT (Fcpos [Work->pivcol] >= 0) ;
    DEBUG6 (("pivcol: "ID" pos "ID" fnr_curr "ID" fncols "ID"\n",
    Work->pivcol, Fcpos [Work->pivcol], Work->fnr_curr, fncols)) ;
    ASSERT (Fcpos [Work->pivcol] <  fncols * Work->fnr_curr) ;
    DEBUG6 (("\npivot col pattern:  fnrows="ID"\n", fnrows)) ;
    for (i = 0 ; i < fnrows ; i++)
    {
    row = Work->Frows [i] ;
    DEBUG6 ((ID" ", row)) ;
    ASSERT (Frpos [row] == i) ;
    ASSERT (NON_PIVOTAL_ROW (row)) ;
    }
    DEBUG6 (("\n")) ;
    ASSERT (Frpos [Work->pivrow] >= 0) ;
    ASSERT (Frpos [Work->pivrow] < fnrows) ;
    ASSERT (Work->Flublock == (Entry *) (Numeric->Memory + E [0])) ;
    ASSERT (Work->Fcblock == Work->Flublock + Work->nb *
    (Work->nb + Work->fnr_curr + Work->fnc_curr)) ;
#endif

    Fcblock = Work->Fcblock ;

    /* ---------------------------------------------------------------------- */
    /* determine the largest actual frontal matrix size (for Info only) */
    /* ---------------------------------------------------------------------- */

    ASSERT (fnrows == Work->fnrows_new + 1) ;
    ASSERT (fncols == Work->fncols_new + 1) ;

    Numeric->maxnrows = MAX (Numeric->maxnrows, fnrows) ;
    Numeric->maxncols = MAX (Numeric->maxncols, fncols) ;

    /* this is safe from integer overflow, since the current frontal matrix
     * is already allocated. */
    Numeric->maxfrsize = MAX (Numeric->maxfrsize, fnrows * fncols) ;

    /* ---------------------------------------------------------------------- */
    /* assemble from prior elements into the current frontal matrix */
    /* ---------------------------------------------------------------------- */

    DEBUG2 (("New assemble start [prior_element:"ID"\n", Work->prior_element)) ;

    /* Currently no rows or columns are marked.  No elements are scanned, */
    /* that is, (ep->next == EMPTY) is true for all elements */

    son_list = 0 ;  /* start creating son_list [ */

    /* ---------------------------------------------------------------------- */
    /* determine if most recent element is Lson or Uson of current front */
    /* ---------------------------------------------------------------------- */

    if (!Work->do_extend)
    {
    prior_Uson = ( Work->pivcol_in_front && !Work->pivrow_in_front) ;
    prior_Lson = (!Work->pivcol_in_front &&  Work->pivrow_in_front) ;
    if (prior_Uson || prior_Lson)
    {
        e = Work->prior_element ;
        if (e != EMPTY)
        {
        ASSERT (E [e]) ;
        p = Memory + E [e] ;
        ep = (Element *) p ;
        ep->next = son_list ;
        son_list = e ;
#ifndef NDEBUG
        DEBUG2 (("e "ID" is Prior son "ID" "ID"\n",
            e, prior_Uson, prior_Lson)) ;
        UMF_dump_element (Numeric, Work, e, FALSE) ;
#endif
        ASSERT (E [e]) ;
        }
    }
    }
    Work->prior_element = EMPTY ;

    /* ---------------------------------------------------------------------- */
    /* SCAN1-row:  scan the element lists of each new row in the pivot col */
    /* and compute the external column degree for each frontal */
    /* ---------------------------------------------------------------------- */

    for (i2 = Work->fscan_row ; i2 < fnrows ; i2++)
    {
    /* Get a row */
    row = Work->NewRows [i2] ;
    if (row < 0) row = FLIP (row) ;
    ASSERT (row >= 0 && row < n_row) ;

    DEBUG6 (("SCAN1-row: "ID"\n", row)) ;
#ifndef NDEBUG
    UMF_dump_rowcol (0, Numeric, Work, row, FALSE) ;
#endif

    ASSERT (NON_PIVOTAL_ROW (row)) ;
    tpi = Row_tuples [row] ;
    if (!tpi) continue ;
    tp = (Tuple *) (Memory + tpi) ;
    tp1 = tp ;
    tp2 = tp ;
    tpend = tp + Row_tlen [row] ;
    for ( ; tp < tpend ; tp++)
    {
        e = tp->e ;
        ASSERT (e > 0 && e <= Work->nel) ;
        if (!E [e]) continue ;  /* element already deallocated */
        f = tp->f ;
        p = Memory + E [e] ;
        ep = (Element *) p ;
        p += UNITS (Element, 1) ;
        Rows = ((Int *) p) + ep->ncols ;
        if (Rows [f] == EMPTY) continue ;   /* row already assembled */
        ASSERT (row == Rows [f]) ;

        if (ep->cdeg < cdeg0)
        {
        /* first time seen in scan1-row */
        ep->cdeg = ep->nrowsleft + cdeg0 ;
        DEBUG6 (("e "ID" First seen: cdeg: "ID" ", e, ep->cdeg-cdeg0)) ;
        ASSERT (ep->ncolsleft > 0 && ep->nrowsleft > 0) ;
        }

        ep->cdeg-- ;    /* decrement external column degree */
        DEBUG6 (("e "ID" New ext col deg: "ID"\n", e, ep->cdeg - cdeg0)) ;

        /* this element is not yet in the new son list */
        if (ep->cdeg == cdeg0 && ep->next == EMPTY)
        {
        /* A new LUson or Uson has been found */
        ep->next = son_list ;
        son_list = e ;
        }

        ASSERT (ep->cdeg >= cdeg0) ;
        *tp2++ = *tp ;  /* leave the tuple in the list */
    }
    Row_tlen [row] = tp2 - tp1 ;
    }

    /* ---------------------------------------------------------------------- */
    /* SCAN1-col:  scan the element lists of each new col in the pivot row */
    /*   and compute the external row degree for each frontal */
    /* ---------------------------------------------------------------------- */

    for (j2 = Work->fscan_col ; j2 < fncols ; j2++)
    {
    /* Get a column */
    col = Work->NewCols [j2] ;
    if (col < 0) col = FLIP (col) ;
    ASSERT (col >= 0 && col < n_col) ;

    DEBUG6 (("SCAN 1-col: "ID"\n", col)) ;
#ifndef NDEBUG
    UMF_dump_rowcol (1, Numeric, Work, col, FALSE) ;
#endif

    ASSERT (NON_PIVOTAL_COL (col)) ;
    tpi = Col_tuples [col] ;
    if (!tpi) continue ;
    tp = (Tuple *) (Memory + tpi) ;
    tp1 = tp ;
    tp2 = tp ;
    tpend = tp + Col_tlen [col] ;
    for ( ; tp < tpend ; tp++)
    {
        e = tp->e ;
        ASSERT (e > 0 && e <= Work->nel) ;
        if (!E [e]) continue ;  /* element already deallocated */
        f = tp->f ;
        p = Memory + E [e] ;
        ep = (Element *) p ;
        p += UNITS (Element, 1) ;
        Cols = (Int *) p ;
        if (Cols [f] == EMPTY) continue ;   /* column already assembled */
        ASSERT (col == Cols [f]) ;

        if (ep->rdeg < rdeg0)
        {
        /* first time seen in scan1-col */
        ep->rdeg = ep->ncolsleft + rdeg0 ;
        DEBUG6 (("e "ID" First seen: rdeg: "ID" ", e, ep->rdeg-rdeg0)) ;
        ASSERT (ep->ncolsleft > 0 && ep->nrowsleft > 0) ;
        }

        ep->rdeg-- ;    /* decrement external row degree */
        DEBUG6 (("e "ID" New ext row degree: "ID"\n", e, ep->rdeg-rdeg0)) ;

        if (ep->rdeg == rdeg0 && ep->next == EMPTY)
        {
        /* A new LUson or Lson has been found */
        ep->next = son_list ;
        son_list = e ;
        }

        ASSERT (ep->rdeg >= rdeg0) ;
        *tp2++ = *tp ;  /* leave the tuple in the list */
    }
    Col_tlen [col] = tp2 - tp1 ;
    }

    /* ---------------------------------------------------------------------- */
    /* assemble new sons via full scans */
    /* ---------------------------------------------------------------------- */

    next = EMPTY ;

    for (e = son_list ; e > 0 ; e = next)
    {
    ASSERT (e > 0 && e <= Work->nel && E [e]) ;
    p = Memory + E [e] ;
    DEBUG2 (("New son: "ID"\n", e)) ;
#ifndef NDEBUG
    UMF_dump_element (Numeric, Work, e, FALSE) ;
#endif
    GET_ELEMENT (ep, p, Cols, Rows, ncols, nrows, S) ;
    nrowsleft = ep->nrowsleft ;
    ncolsleft = ep->ncolsleft ;
    next = ep->next ;
    ep->next = EMPTY ;

    extrdeg = (ep->rdeg < rdeg0) ? ncolsleft : (ep->rdeg - rdeg0) ;
    extcdeg = (ep->cdeg < cdeg0) ? nrowsleft : (ep->cdeg - cdeg0) ;
    ncols_to_assemble = ncolsleft - extrdeg ;
    nrows_to_assemble = nrowsleft - extcdeg ;
    DEBUG2 (("extrdeg "ID" extcdeg "ID"\n", extrdeg, extcdeg)) ;

    if (extrdeg == 0 && extcdeg == 0)
    {

        /* -------------------------------------------------------------- */
        /* this is an LUson - assemble an entire contribution block */
        /* -------------------------------------------------------------- */

        DEBUG6 (("LUson found: "ID"\n", e)) ;

        if (nrows == nrowsleft)
        {
        /* ---------------------------------------------------------- */
        /* no rows assembled out of this LUson yet */
        /* ---------------------------------------------------------- */

        /* compute the compressed column offset vector*/
        /* [ use Wm [0..nrows-1] for offsets */
#pragma ivdep
        for (i = 0 ; i < nrows ; i++)
        {
            row = Rows [i] ;
            Row_degree [row] -= ncolsleft ;
            Wm [i] = Frpos [row] ;
        }

        if (ncols == ncolsleft)
        {
            /* ------------------------------------------------------ */
            /* no rows or cols assembled out of LUson yet */
            /* ------------------------------------------------------ */

            for (j = 0 ; j < ncols ; j++)
            {
            col = Cols [j] ;
#ifndef FIXQ
            Col_degree [col] -= nrowsleft ;
#endif
            Fcol = Fcblock + Fcpos [col] ;
#pragma ivdep
            for (i = 0 ; i < nrows ; i++)
            {
                /* Fcol [Wm [i]] += S [i] ; */
                ASSEMBLE (Fcol [Wm [i]], S [i]) ;
            }
            S += nrows ;
            }


        }
        else
        {
            /* ------------------------------------------------------ */
            /* only cols have been assembled out of LUson */
            /* ------------------------------------------------------ */

            for (j = 0 ; j < ncols ; j++)
            {
            col = Cols [j] ;
            if (col >= 0)
            {
#ifndef FIXQ
                Col_degree [col] -= nrowsleft ;
#endif
                Fcol = Fcblock + Fcpos [col] ;
#pragma ivdep
                for (i = 0 ; i < nrows ; i++)
                {
                /* Fcol [Wm [i]] += S [i] ; */
                ASSEMBLE (Fcol [Wm [i]], S [i]) ;
                }
            }
            S += nrows ;
            }

        }
        /* ] done using Wm [0..nrows-1] for offsets */
        }
        else
        {
        /* ---------------------------------------------------------- */
        /* some rows have been assembled out of this LUson */
        /* ---------------------------------------------------------- */

        /* compute the compressed column offset vector*/
        /* [ use Woo,Wm [0..nrowsleft-1] for offsets */
        ngetrows = 0 ;
        for (i = 0 ; i < nrows ; i++)
        {
            row = Rows [i] ;
            if (row >= 0)
            {
            Row_degree [row] -= ncolsleft ;
            Woo [ngetrows] = i ;
            Wm [ngetrows++] = Frpos [row] ;
            }
        }
        ASSERT (ngetrows == nrowsleft) ;

        if (ncols == ncolsleft)
        {
            /* ------------------------------------------------------ */
            /* only rows have been assembled out of this LUson */
            /* ------------------------------------------------------ */

            for (j = 0 ; j < ncols ; j++)
            {
            col = Cols [j] ;
#ifndef FIXQ
            Col_degree [col] -= nrowsleft ;
#endif
            Fcol = Fcblock + Fcpos [col] ;
#pragma ivdep
            for (i = 0 ; i < nrowsleft ; i++)
            {
                /* Fcol [Wm [i]] += S [Woo [i]] ; */
                ASSEMBLE (Fcol [Wm [i]], S [Woo [i]]) ;
            }
            S += nrows ;
            }

        }
        else
        {
            /* ------------------------------------------------------ */
            /* both rows and columns have been assembled out of LUson */
            /* ------------------------------------------------------ */

            for (j = 0 ; j < ncols ; j++)
            {
            col = Cols [j] ;
            if (col >= 0)
            {
#ifndef FIXQ
                Col_degree [col] -= nrowsleft ;
#endif
                Fcol = Fcblock + Fcpos [col] ;
#pragma ivdep
                for (i = 0 ; i < nrowsleft ; i++)
                {
                /* Fcol [Wm [i]] += S [Woo [i]] ; */
                ASSEMBLE (Fcol [Wm [i]], S [Woo [i]]) ;
                }
            }
            S += nrows ;
            }

        }
        /* ] done using Woo,Wm [0..nrowsleft-1] */
        }

        /* deallocate the element: remove from ordered list */
        UMF_mem_free_tail_block (Numeric, E [e]) ;
        E [e] = 0 ;

    }
    else if (extcdeg == 0)
    {

        /* -------------------------------------------------------------- */
        /* this is a Uson - assemble all possible columns */
        /* -------------------------------------------------------------- */

        DEBUG6 (("New USON: "ID"\n", e)) ;
        ASSERT (extrdeg > 0) ;

        DEBUG6 (("New uson "ID" cols to do "ID"\n", e, ncols_to_assemble)) ;

        if (ncols_to_assemble > 0)
        {

        Int skip = FALSE ;
        if (ncols_to_assemble * 16 < ncols && nrows == 1)
        {
            /* this is a tall and thin frontal matrix consisting of
             * only one column (most likely an original column). Do
             * not assemble it.   It cannot be the pivot column, since
             * the pivot column element would be an LU son, not an Lson,
             * of the current frontal matrix. */
            ASSERT (nrowsleft == 1) ;
            ASSERT (Rows [0] >= 0 && Rows [0] < Work->n_row) ;
            skip = TRUE ;
            Work->any_skip = TRUE ;
        }

        if (!skip)
        {

            if (nrows == nrowsleft)
            {
            /* -------------------------------------------------- */
            /* no rows have been assembled out of this Uson yet */
            /* -------------------------------------------------- */

            /* compute the compressed column offset vector */
            /* [ use Wm [0..nrows-1] for offsets */
#pragma ivdep
            for (i = 0 ; i < nrows ; i++)
            {
                row = Rows [i] ;
                ASSERT (row >= 0 && row < n_row) ;
                Row_degree [row] -= ncols_to_assemble ;
                Wm [i] = Frpos [row] ;
            }

            for (j = 0 ; j < ncols ; j++)
            {
                col = Cols [j] ;
                if ((col >= 0) && (Fcpos [col] >= 0))
                {
#ifndef FIXQ
                Col_degree [col] -= nrowsleft ;
#endif
                Fcol = Fcblock + Fcpos [col] ;
#pragma ivdep
                for (i = 0 ; i < nrows ; i++)
                {
                    /* Fcol [Wm [i]] += S [i] ; */
                    ASSEMBLE (Fcol [Wm [i]], S [i]) ;
                }
                /* flag the column as assembled from Uson */
                Cols [j] = EMPTY ;
                }
                S += nrows ;
            }


            /* ] done using Wm [0..nrows-1] for offsets */
            }
            else
            {
            /* -------------------------------------------------- */
            /* some rows have been assembled out of this Uson */
            /* -------------------------------------------------- */

            /* compute the compressed column offset vector*/
            /* [ use Woo,Wm [0..nrows-1] for offsets */
            ngetrows = 0 ;
            for (i = 0 ; i < nrows ; i++)
            {
                row = Rows [i] ;
                if (row >= 0)
                {
                Row_degree [row] -= ncols_to_assemble ;
                ASSERT (row < n_row && Frpos [row] >= 0) ;
                Woo [ngetrows] = i ;
                Wm [ngetrows++] = Frpos [row] ;
                }
            }
            ASSERT (ngetrows == nrowsleft) ;

            for (j = 0 ; j < ncols ; j++)
            {
                col = Cols [j] ;
                if ((col >= 0) && (Fcpos [col] >= 0))
                {
#ifndef FIXQ
                Col_degree [col] -= nrowsleft ;
#endif
                Fcol = Fcblock + Fcpos [col] ;
#pragma ivdep
                for (i = 0 ; i < nrowsleft ; i++)
                {
                    /* Fcol [Wm [i]] += S [Woo [i]] ; */
                    ASSEMBLE (Fcol [Wm [i]], S [Woo [i]]) ;
                }
                /* flag the column as assembled from Uson */
                Cols [j] = EMPTY ;
                }
                S += nrows ;
            }

            /* ] done using Woo,Wm */
            }
            ep->ncolsleft = extrdeg ;
        }
        }

    }
    else
    {

        /* -------------------------------------------------------------- */
        /* this is an Lson - assemble all possible rows */
        /* -------------------------------------------------------------- */

        DEBUG6 (("New LSON: "ID"\n", e)) ;
        ASSERT (extrdeg == 0 && extcdeg > 0) ;

        DEBUG6 (("New lson "ID" rows to do "ID"\n", e, nrows_to_assemble)) ;

        if (nrows_to_assemble > 0)
        {

        Int skip = FALSE ;
        if (nrows_to_assemble * 16 < nrows && ncols == 1)
        {
            /* this is a tall and thin frontal matrix consisting of
             * only one column (most likely an original column). Do
             * not assemble it.   It cannot be the pivot column, since
             * the pivot column element would be an LU son, not an Lson,
             * of the current frontal matrix. */
            ASSERT (ncolsleft == 1) ;
            ASSERT (Cols [0] >= 0 && Cols [0] < Work->n_col) ;
            Work->any_skip = TRUE ;
            skip = TRUE ;
        }

        if (!skip)
        {

            /* compute the compressed column offset vector */
            /* [ use Woo,Wm [0..nrows-1] for offsets */
            ngetrows = 0 ;
            for (i = 0 ; i < nrows ; i++)
            {
            row = Rows [i] ;
            if ((row >= 0) && (Frpos [row] >= 0))
            {
                ASSERT (row < n_row) ;
                Row_degree [row] -= ncolsleft ;
                Woo [ngetrows] = i ;
                Wm [ngetrows++] = Frpos [row] ;
                /* flag the row as assembled from the Lson */
                Rows [i] = EMPTY ;
            }
            }
            ASSERT (nrowsleft - ngetrows == extcdeg) ;
            ASSERT (ngetrows == nrows_to_assemble) ;

            if (ncols == ncolsleft)
            {
            /* -------------------------------------------------- */
            /* no columns assembled out this Lson yet */
            /* -------------------------------------------------- */

            for (j = 0 ; j < ncols ; j++)
            {
                col = Cols [j] ;
                ASSERT (col >= 0 && col < n_col) ;
#ifndef FIXQ
                Col_degree [col] -= nrows_to_assemble ;
#endif
                Fcol = Fcblock + Fcpos [col] ;
#pragma ivdep
                for (i = 0 ; i < nrows_to_assemble ; i++)
                {
                /* Fcol [Wm [i]] += S [Woo [i]] ; */
                ASSEMBLE (Fcol [Wm [i]], S [Woo [i]]) ;
                }
                S += nrows ;
            }


            }
            else
            {
            /* -------------------------------------------------- */
            /* some columns have been assembled out of this Lson */
            /* -------------------------------------------------- */

            for (j = 0 ; j < ncols ; j++)
            {
                col = Cols [j] ;
                ASSERT (col < n_col) ;
                if (col >= 0)
                {
#ifndef FIXQ
                Col_degree [col] -= nrows_to_assemble ;
#endif
                Fcol = Fcblock + Fcpos [col] ;
#pragma ivdep
                for (i = 0 ; i < nrows_to_assemble ; i++)
                {
                    /* Fcol [Wm [i]] += S [Woo [i]] ; */
                    ASSEMBLE (Fcol [Wm [i]], S [Woo [i]]) ;
                }
                }
                S += nrows ;
            }

            }

            /* ] done using Woo,Wm */

            ep->nrowsleft = extcdeg ;
        }
        }
    }
    }

    /* Note that garbage collection, and build tuples */
    /* both destroy the son list. */

    /* ] son_list now empty */

    /* ---------------------------------------------------------------------- */
    /* If frontal matrix extended, assemble old L/Usons from new rows/cols */
    /* ---------------------------------------------------------------------- */

    /* ---------------------------------------------------------------------- */
    /* SCAN2-row:  assemble rows of old Lsons from the new rows */
    /* ---------------------------------------------------------------------- */

#ifndef NDEBUG
    DEBUG7 (("Current frontal matrix: (prior to scan2-row)\n")) ;
    UMF_dump_current_front (Numeric, Work, TRUE) ;
#endif

    /* rescan the pivot row */
    if (Work->any_skip)
    {
    row_assemble (Work->pivrow, Numeric, Work) ;
    }

    if (Work->do_scan2row)
    {
    for (i2 = Work->fscan_row ; i2 < fnrows ; i2++)
    {
        /* Get a row */
        row = Work->NewRows [i2] ;
        if (row < 0) row = FLIP (row) ;
        ASSERT (row >= 0 && row < n_row) ;
        if (!(row == Work->pivrow && Work->any_skip))
        {
        /* assemble it */
        row_assemble (row, Numeric, Work) ;
        }
    }
    }

    /* ---------------------------------------------------------------------- */
    /* SCAN2-col:  assemble columns of old Usons from the new columns */
    /* ---------------------------------------------------------------------- */

#ifndef NDEBUG
    DEBUG7 (("Current frontal matrix: (prior to scan2-col)\n")) ;
    UMF_dump_current_front (Numeric, Work, TRUE) ;
#endif

    /* rescan the pivot col */
    if (Work->any_skip)
    {
    col_assemble (Work->pivcol, Numeric, Work) ;
    }

    if (Work->do_scan2col)
    {

    for (j2 = Work->fscan_col ; j2 < fncols ; j2++)
    {
        /* Get a column */
        col = Work->NewCols [j2] ;
        if (col < 0) col = FLIP (col) ;
        ASSERT (col >= 0 && col < n_col) ;
        if (!(col == Work->pivcol && Work->any_skip))
        {
        /* assemble it */
        col_assemble (col, Numeric, Work) ;
        }
    }
    }

    /* ---------------------------------------------------------------------- */
    /* done.  the remainder of this routine is used only when in debug mode */
    /* ---------------------------------------------------------------------- */

#ifndef NDEBUG

    /* ---------------------------------------------------------------------- */
    /* when debugging: make sure the assembly did everything that it could */
    /* ---------------------------------------------------------------------- */

    DEBUG3 (("::Assemble done\n")) ;

    for (i2 = 0 ; i2 < fnrows ; i2++)
    {
    /* Get a row */
    row = Work->Frows [i2] ;
    ASSERT (row >= 0 && row < n_row) ;

    DEBUG6 (("DEBUG SCAN 1: "ID"\n", row)) ;
    UMF_dump_rowcol (0, Numeric, Work, row, TRUE) ;

    ASSERT (NON_PIVOTAL_ROW (row)) ;
    tpi = Row_tuples [row] ;
    if (!tpi) continue ;
    tp = (Tuple *) (Memory + tpi) ;
    tpend = tp + Row_tlen [row] ;
    for ( ; tp < tpend ; tp++)
    {
        e = tp->e ;
        ASSERT (e > 0 && e <= Work->nel) ;
        if (!E [e]) continue ;  /* element already deallocated */
        f = tp->f ;
        p = Memory + E [e] ;
        ep = (Element *) p ;
        p += UNITS (Element, 1) ;
        Cols = (Int *) p ;
        Rows = ((Int *) p) + ep->ncols ;
        if (Rows [f] == EMPTY) continue ;   /* row already assembled */
        ASSERT (row == Rows [f]) ;
        extrdeg = (ep->rdeg < rdeg0) ? ep->ncolsleft : (ep->rdeg - rdeg0) ;
        extcdeg = (ep->cdeg < cdeg0) ? ep->nrowsleft : (ep->cdeg - cdeg0) ;
        DEBUG6 ((
        "e "ID" After assembly ext row deg: "ID" ext col degree "ID"\n",
        e, extrdeg, extcdeg)) ;

        if (Work->any_skip)
        {
        /* no Lsons in any row, except for very tall and thin ones */
        ASSERT (extrdeg >= 0) ;
        if (extrdeg == 0)
        {
            /* this is an unassemble Lson */
            ASSERT (ep->ncols == 1) ;
            ASSERT (ep->ncolsleft == 1) ;
            col = Cols [0] ;
            ASSERT (col != Work->pivcol) ;
        }
        }
        else
        {
        /* no Lsons in any row */
        ASSERT (extrdeg > 0) ;
        /* Uson external row degree is = number of cols left */
        ASSERT (IMPLIES (extcdeg == 0, extrdeg == ep->ncolsleft)) ;
        }
    }
    }

    /* ---------------------------------------------------------------------- */

    for (j2 = 0 ; j2 < fncols ; j2++)
    {
    /* Get a column */
    col = Work->Fcols [j2] ;
    ASSERT (col >= 0 && col < n_col) ;

    DEBUG6 (("DEBUG SCAN 2: "ID"\n", col)) ;
#ifndef FIXQ
    UMF_dump_rowcol (1, Numeric, Work, col, TRUE) ;
#else
    UMF_dump_rowcol (1, Numeric, Work, col, FALSE) ;
#endif

    ASSERT (NON_PIVOTAL_COL (col)) ;
    tpi = Col_tuples [col] ;
    if (!tpi) continue ;
    tp = (Tuple *) (Memory + tpi) ;
    tpend = tp + Col_tlen [col] ;
    for ( ; tp < tpend ; tp++)
    {
        e = tp->e ;
        ASSERT (e > 0 && e <= Work->nel) ;
        if (!E [e]) continue ;  /* element already deallocated */
        f = tp->f ;
        p = Memory + E [e] ;
        ep = (Element *) p ;
        p += UNITS (Element, 1) ;
        Cols = (Int *) p ;
        Rows = ((Int *) p) + ep->ncols ;
        if (Cols [f] == EMPTY) continue ;   /* column already assembled */
        ASSERT (col == Cols [f]) ;
        extrdeg = (ep->rdeg < rdeg0) ? ep->ncolsleft : (ep->rdeg - rdeg0) ;
        extcdeg = (ep->cdeg < cdeg0) ? ep->nrowsleft : (ep->cdeg - cdeg0) ;
        DEBUG6 (("e "ID" After assembly ext col deg: "ID"\n", e, extcdeg)) ;

        if (Work->any_skip)
        {
        /* no Usons in any column, except for very tall and thin ones */
        ASSERT (extcdeg >= 0) ;
        if (extcdeg == 0)
        {
            /* this is an unassemble Uson */
            ASSERT (ep->nrows == 1) ;
            ASSERT (ep->nrowsleft == 1) ;
            row = Rows [0] ;
            ASSERT (row != Work->pivrow) ;
        }
        }
        else
        {
        /* no Usons in any column */
        ASSERT (extcdeg > 0) ;
        /* Lson external column degree is = number of rows left */
        ASSERT (IMPLIES (extrdeg == 0, extcdeg == ep->nrowsleft)) ;
        }
    }
    }
#endif /* NDEBUG */
}