The R Project SVN R-packages

/* ========================================================================== */
/* === UMF_row_search ======================================================= */
/* ========================================================================== */

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

/*
    Find two candidate pivot rows in a column: the best one in the front,
    and the best one not in the front.  Return the two pivot row patterns and
    their exact degrees.  Called by UMF_local_search.

    Returns UMFPACK_OK if successful, or UMFPACK_WARNING_singular_matrix or
    UMFPACK_ERROR_different_pattern if not.

*/

#include "umf_internal.h"
#include "umf_row_search.h"

GLOBAL Int UMF_row_search
(
    NumericType *Numeric,
    WorkType *Work,
    SymbolicType *Symbolic,
    Int cdeg0,          /* length of column in Front */
    Int cdeg1,          /* length of column outside Front */
    const Int Pattern [ ],  /* pattern of column, Pattern [0..cdeg1 -1] */
    const Int Pos [ ],      /* Pos [Pattern [0..cdeg1 -1]] = 0..cdeg1 -1 */
    Int pivrow [2],     /* pivrow [IN] and pivrow [OUT] */
    Int rdeg [2],       /* rdeg [IN] and rdeg [OUT] */
    Int W_i [ ],        /* pattern of pivrow [IN], */
                /* either Fcols or Woi */
    Int W_o [ ],        /* pattern of pivrow [OUT], */
                /* either Wio or Woo */
    Int prior_pivrow [2],   /* the two other rows just scanned, if any */
    const Entry Wxy [ ],    /* numerical values Wxy [0..cdeg1-1],
                   either Wx or Wy */

    Int pivcol,         /* the candidate column being searched */
    Int freebie [ ]
)
{

    /* ---------------------------------------------------------------------- */
    /* local variables */
    /* ---------------------------------------------------------------------- */

    double maxval, toler, toler2, value, pivot [2] ;
    Int i, row, deg, col, *Frpos, fnrows, *E, j, ncols, *Cols, *Rows,
    e, f, Wrpflag, *Fcpos, fncols, tpi, max_rdeg, nans_in_col, was_offdiag,
    diag_row, prefer_diagonal, *Wrp, found, *Diagonal_map ;
    Tuple *tp, *tpend, *tp1, *tp2 ;
    Unit *Memory, *p ;
    Element *ep ;
    Int *Row_tuples, *Row_degree, *Row_tlen ;

#ifndef NDEBUG
    Int *Col_degree ;
    DEBUG2 (("Row_search:\n")) ;
    for (i = 0 ; i < cdeg1 ; i++)
    {
    row = Pattern [i] ;
    DEBUG4 (("   row: "ID"\n", row)) ;
    ASSERT (row >= 0 && row < Numeric->n_row) ;
    ASSERT (i == Pos [row]) ;
    }
    /* If row is not in Pattern [0..cdeg1-1], then Pos [row] == EMPTY */
    if (UMF_debug > 0 || Numeric->n_row < 1000)
    {
    Int cnt = cdeg1 ;
    DEBUG4 (("Scan all rows:\n")) ;
    for (row = 0 ; row < Numeric->n_row ; row++)
    {
        if (Pos [row] < 0)
        {
        cnt++ ;
        }
        else
        {
        DEBUG4 (("   row: "ID" pos "ID"\n", row, Pos [row])) ;
        }
    }
    ASSERT (cnt == Numeric->n_row) ;
    }
    Col_degree = Numeric->Cperm ;   /* for NON_PIVOTAL_COL macro only */
    ASSERT (pivcol >= 0 && pivcol < Work->n_col) ;
    ASSERT (NON_PIVOTAL_COL (pivcol)) ;
#endif

    pivot [IN] = 0. ;
    pivot [OUT] = 0. ;

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

    Row_degree = Numeric->Rperm ;
    Row_tuples = Numeric->Uip ;
    Row_tlen   = Numeric->Uilen ;
    Wrp = Work->Wrp ;
    Frpos = Work->Frpos ;
    E = Work->E ;
    Memory = Numeric->Memory ;
    fnrows = Work->fnrows ;

    prefer_diagonal = Symbolic->prefer_diagonal ;
    Diagonal_map = Work->Diagonal_map ;

    if (Diagonal_map)
    {
    diag_row = Diagonal_map [pivcol] ;
    was_offdiag = diag_row < 0 ;
    if (was_offdiag)
    {
        /* the "diagonal" entry in this column was permuted here by an
         * earlier pivot choice.  The tighter off-diagonal tolerance will
         * be used instead of the symmetric tolerance. */
        diag_row = FLIP (diag_row) ;
    }
    ASSERT (diag_row >= 0 && diag_row < Numeric->n_row) ;
    }
    else
    {
    diag_row = EMPTY ;  /* unused */
    was_offdiag = EMPTY ;   /* unused */
    }

    /* pivot row degree cannot exceed max_rdeg */
    max_rdeg = Work->fncols_max ;

    /* ---------------------------------------------------------------------- */
    /* scan pivot column for candidate rows */
    /* ---------------------------------------------------------------------- */

    maxval = 0.0 ;
    nans_in_col = FALSE ;

    for (i = 0 ; i < cdeg1 ; i++)
    {
    APPROX_ABS (value, Wxy [i]) ;
    if (SCALAR_IS_NAN (value))
    {
        nans_in_col = TRUE ;
        maxval = value ;
        break ;
    }
    /* This test can now ignore the NaN case: */
    maxval = MAX (maxval, value) ;
    }

    /* if maxval is zero, the matrix is numerically singular */

    toler = Numeric->relpt * maxval ;
    toler2 = Numeric->relpt2 * maxval ;
    toler2 = was_offdiag ? toler : toler2 ;

    DEBUG5 (("Row_search begins [ maxval %g toler %g %g\n",
    maxval, toler, toler2)) ;
    if (SCALAR_IS_NAN (toler) || SCALAR_IS_NAN (toler2))
    {
    nans_in_col = TRUE ;
    }

    if (!nans_in_col)
    {

    /* look for the diagonal entry, if it exists */
    found = FALSE ;
    ASSERT (!SCALAR_IS_NAN (toler)) ;

    if (prefer_diagonal)
    {
        ASSERT (diag_row != EMPTY) ;
        i = Pos [diag_row] ;
        if (i >= 0)
        {
        double a ;
        ASSERT (i < cdeg1) ;
        ASSERT (diag_row == Pattern [i]) ;

        APPROX_ABS (a, Wxy [i]) ;

        ASSERT (!SCALAR_IS_NAN (a)) ;
        ASSERT (!SCALAR_IS_NAN (toler2)) ;

        if (SCALAR_IS_NONZERO (a) && a >= toler2)
        {
            /* found it! */
            DEBUG3 (("Symmetric pivot: "ID" "ID"\n", pivcol, diag_row));
            found = TRUE ;
            if (Frpos [diag_row] >= 0 && Frpos [diag_row] < fnrows)
            {
            pivrow [IN] = diag_row ;
            pivrow [OUT] = EMPTY ;
            }
            else
            {
            pivrow [IN] = EMPTY ;
            pivrow [OUT] = diag_row ;
            }
        }
        }
    }

    /* either no diagonal found, or we didn't look for it */
    if (!found)
    {
        if (cdeg0 > 0)
        {

        /* this is a column in the front */
        for (i = 0 ; i < cdeg0 ; i++)
        {
            double a ;
            APPROX_ABS (a, Wxy [i]) ;
            ASSERT (!SCALAR_IS_NAN (a)) ;
            ASSERT (!SCALAR_IS_NAN (toler)) ;
            if (SCALAR_IS_NONZERO (a) && a >= toler)
            {
            row = Pattern [i] ;
            deg = Row_degree [row] ;
#ifndef NDEBUG
            DEBUG6 ((ID" candidate row "ID" deg "ID" absval %g\n",
                i, row, deg, a)) ;
            UMF_dump_rowcol (0, Numeric, Work, row, TRUE) ;
#endif
            ASSERT (Frpos [row] >= 0 && Frpos [row] < fnrows) ;
            ASSERT (Frpos [row] == i) ;
            /* row is in the current front */
            DEBUG4 ((" in front\n")) ;
            if (deg < rdeg [IN]
                /* break ties by picking the largest entry: */
                   || (deg == rdeg [IN] && a > pivot [IN])
                /* break ties by picking the diagonal entry: */
                /* || (deg == rdeg [IN] && row == diag_row) */
               )
            {
                /* best row in front, so far */
                pivrow [IN] = row ;
                rdeg [IN] = deg ;
                pivot [IN] = a ;
            }
            }
        }
        for ( ; i < cdeg1 ; i++)
        {
            double a ;
            APPROX_ABS (a, Wxy [i]) ;
            ASSERT (!SCALAR_IS_NAN (a)) ;
            ASSERT (!SCALAR_IS_NAN (toler)) ;
            if (SCALAR_IS_NONZERO (a) && a >= toler)
            {
            row = Pattern [i] ;
            deg = Row_degree [row] ;
#ifndef NDEBUG
            DEBUG6 ((ID" candidate row "ID" deg "ID" absval %g\n",
                i, row, deg, a)) ;
            UMF_dump_rowcol (0, Numeric, Work, row, TRUE) ;
#endif
            ASSERT (Frpos [row] == i) ;
            /* row is not in the current front */
            DEBUG4 ((" NOT in front\n")) ;
            if (deg < rdeg [OUT]
                /* break ties by picking the largest entry: */
                   || (deg == rdeg [OUT] && a > pivot [OUT])
                /* break ties by picking the diagonal entry: */
                /* || (deg == rdeg [OUT] && row == diag_row) */
               )
            {
                /* best row not in front, so far */
                pivrow [OUT] = row ;
                rdeg [OUT] = deg ;
                pivot [OUT] = a ;
            }
            }
        }

        }
        else
        {

        /* this column is not in the front */
        for (i = 0 ; i < cdeg1 ; i++)
        {
            double a ;
            APPROX_ABS (a, Wxy [i]) ;
            ASSERT (!SCALAR_IS_NAN (a)) ;
            ASSERT (!SCALAR_IS_NAN (toler)) ;
            if (SCALAR_IS_NONZERO (a) && a >= toler)
            {
            row = Pattern [i] ;
            deg = Row_degree [row] ;
#ifndef NDEBUG
            DEBUG6 ((ID" candidate row "ID" deg "ID" absval %g\n",
                i, row, deg, a)) ;
            UMF_dump_rowcol (0, Numeric, Work, row, TRUE) ;
#endif
            if (Frpos [row] >= 0 && Frpos [row] < fnrows)
            {
                /* row is in the current front */
                DEBUG4 ((" in front\n")) ;
                if (deg < rdeg [IN]
                /* break ties by picking the largest entry: */
                   || (deg == rdeg [IN] && a > pivot [IN])
                /* break ties by picking the diagonal entry: */
                /* || (deg == rdeg [IN] && row == diag_row) */
                   )
                {
                /* best row in front, so far */
                pivrow [IN] = row ;
                rdeg [IN] = deg ;
                pivot [IN] = a ;
                }
            }
            else
            {
                /* row is not in the current front */
                DEBUG4 ((" NOT in front\n")) ;
                if (deg < rdeg [OUT]
                /* break ties by picking the largest entry: */
                   || (deg == rdeg[OUT] && a > pivot [OUT])
                /* break ties by picking the diagonal entry: */
                /* || (deg == rdeg[OUT] && row == diag_row) */
                   )
                {
                /* best row not in front, so far */
                pivrow [OUT] = row ;
                rdeg [OUT] = deg ;
                pivot [OUT] = a ;
                }
            }
            }
        }
        }
    }
    }

    /* ---------------------------------------------------------------------- */
    /* NaN handling */
    /* ---------------------------------------------------------------------- */

    /* if cdeg1 > 0 then we must have found a pivot row ... unless NaN's */
    /* exist.  Try with no numerical tests if no pivot found. */

    if (cdeg1 > 0 && pivrow [IN] == EMPTY && pivrow [OUT] == EMPTY)
    {
    /* cleanup for the NaN case */
    DEBUG0 (("Found a NaN in pivot column!\n")) ;

    /* grab the first entry in the pivot column, ignoring degree, */
    /* numerical stability, and symmetric preference */
    row = Pattern [0] ;
    deg = Row_degree [row] ;
    if (Frpos [row] >= 0 && Frpos [row] < fnrows)
    {
        /* row is in the current front */
        DEBUG4 ((" in front\n")) ;
        pivrow [IN] = row ;
        rdeg [IN] = deg ;
    }
    else
    {
        /* row is not in the current front */
        DEBUG4 ((" NOT in front\n")) ;
        pivrow [OUT] = row ;
        rdeg [OUT] = deg ;
    }

    /* We are now guaranteed to have a pivot, no matter how broken */
    /* (non-IEEE compliant) the underlying numerical operators are. */
    /* This is particularly a problem for Microsoft compilers (they do */
    /* not handle NaN's properly). Now try to find a sparser pivot, if */
    /* possible. */

    for (i = 1 ; i < cdeg1 ; i++)
    {
        row = Pattern [i] ;
        deg = Row_degree [row] ;

        if (Frpos [row] >= 0 && Frpos [row] < fnrows)
        {
        /* row is in the current front */
        DEBUG4 ((" in front\n")) ;
        if (deg < rdeg [IN] || (deg == rdeg [IN] && row == diag_row))
        {
            /* best row in front, so far */
            pivrow [IN] = row ;
            rdeg [IN] = deg ;
        }
        }
        else
        {
        /* row is not in the current front */
        DEBUG4 ((" NOT in front\n")) ;
        if (deg < rdeg [OUT] || (deg == rdeg [OUT] && row == diag_row))
        {
            /* best row not in front, so far */
            pivrow [OUT] = row ;
            rdeg [OUT] = deg ;
        }
        }
    }
    }

    /* We found a pivot if there are entries (even zero ones) in pivot col */
    ASSERT (IMPLIES (cdeg1 > 0, pivrow[IN] != EMPTY || pivrow[OUT] != EMPTY)) ;

    /* If there are no entries in the pivot column, then no pivot is found */
    ASSERT (IMPLIES (cdeg1 == 0, pivrow[IN] == EMPTY && pivrow[OUT] == EMPTY)) ;

    /* ---------------------------------------------------------------------- */
    /* check for singular matrix */
    /* ---------------------------------------------------------------------- */

    if (cdeg1  == 0)
    {
    if (fnrows > 0)
    {
        /*
        Get the pivrow [OUT][IN] from the current front.
        The frontal matrix looks like this:

            pivcol[OUT]
            |
            v
        x x x x 0   <- so grab this row as the pivrow [OUT][IN].
        x x x x 0
        x x x x 0
        0 0 0 0 0

        The current frontal matrix has some rows in it.  The degree
        of the pivcol[OUT] is zero.  The column is empty, and the
        current front does not contribute to it.

        */
        pivrow [IN] = Work->Frows [0] ;
        DEBUGm4 (("Got zero pivrow[OUT][IN] "ID" from current front\n",
        pivrow [IN])) ;
    }
    else
    {

        /*
        Get a pivot row from the row-merge tree, use as
        pivrow [OUT][OUT].   pivrow [IN] remains EMPTY.
        This can only happen if the current front is 0-by-0.
        */

        Int *Front_leftmostdesc, *Front_1strow, *Front_new1strow, row1,
        row2, fleftmost, nfr, n_row, frontid ;

        ASSERT (Work->fncols == 0) ;

        Front_leftmostdesc = Symbolic->Front_leftmostdesc ;
        Front_1strow = Symbolic->Front_1strow ;
        Front_new1strow = Work->Front_new1strow ;
        nfr = Symbolic->nfr ;
        n_row = Numeric->n_row ;
        frontid = Work->frontid ;

        DEBUGm4 (("Note: pivcol: "ID" is empty front "ID"\n",
        pivcol, frontid)) ;
#ifndef NDEBUG
        DEBUG1 (("Calling dump rowmerge\n")) ;
        UMF_dump_rowmerge (Numeric, Symbolic, Work) ;
#endif

        /* Row-merge set is the non-pivotal rows in the range */
        /* Front_new1strow [Front_leftmostdesc [frontid]] to */
        /* Front_1strow [frontid+1] - 1. */
        /* If this is empty, then use the empty rows, in the range */
        /* Front_new1strow [nfr] to n_row-1. */
        /* If this too is empty, then pivrow [OUT] will be empty. */
        /* In both cases, update Front_new1strow [...]. */

        fleftmost = Front_leftmostdesc [frontid] ;
        row1 = Front_new1strow [fleftmost] ;
        row2 = Front_1strow [frontid+1] - 1 ;
        DEBUG1 (("Leftmost: "ID" Rows ["ID" to "ID"] srch ["ID" to "ID"]\n",
        fleftmost, Front_1strow [frontid], row2, row1, row2)) ;

        /* look in the range row1 ... row2 */
        for (row = row1 ; row <= row2 ; row++)
        {
        DEBUG3 (("   Row: "ID"\n", row)) ;
        if (NON_PIVOTAL_ROW (row))
        {
            /* found it */
            DEBUG3 (("   Row: "ID" found\n", row)) ;
            ASSERT (Frpos [row] == EMPTY) ;
            pivrow [OUT] = row ;
            DEBUGm4 (("got row merge pivrow %d\n", pivrow [OUT])) ;
            break ;
        }
        }
        Front_new1strow [fleftmost] = row ;

        if (pivrow [OUT] == EMPTY)
        {
        /* not found, look in empty row set in "dummy" front */
        row1 = Front_new1strow [nfr] ;
        row2 = n_row-1 ;
        DEBUG3 (("Empty: "ID" Rows ["ID" to "ID"] srch["ID" to "ID"]\n",
            nfr, Front_1strow [nfr], row2, row1, row2)) ;

        /* look in the range row1 ... row2 */
        for (row = row1 ; row <= row2 ; row++)
        {
            DEBUG3 (("   Empty Row: "ID"\n", row)) ;
            if (NON_PIVOTAL_ROW (row))
            {
            /* found it */
            DEBUG3 (("   Empty Row: "ID" found\n", row)) ;
            ASSERT (Frpos [row] == EMPTY) ;
            pivrow [OUT] = row ;
            DEBUGm4 (("got dummy row pivrow %d\n", pivrow [OUT])) ;
            break ;
            }
        }
        Front_new1strow [nfr] = row ;
        }

        if (pivrow [OUT] == EMPTY)
        {
        /* Row-merge set is empty.  We can just discard */
        /* the candidate pivot column. */
        DEBUG0 (("Note: row-merge set empty\n")) ;
        DEBUGm4 (("got no pivrow \n")) ;
        return (UMFPACK_WARNING_singular_matrix) ;
        }
    }
    }

    /* ---------------------------------------------------------------------- */
    /* construct the candidate row in the front, if any */
    /* ---------------------------------------------------------------------- */

#ifndef NDEBUG
    /* check Wrp */
    ASSERT (Work->Wrpflag > 0) ;
    if (UMF_debug > 0 || Work->n_col < 1000)
    {
    for (i = 0 ; i < Work->n_col ; i++)
    {
        ASSERT (Wrp [i] < Work->Wrpflag) ;
    }
    }
#endif

#ifndef NDEBUG
    DEBUG4 (("pivrow [IN]: "ID"\n", pivrow [IN])) ;
    UMF_dump_rowcol (0, Numeric, Work, pivrow [IN], TRUE) ;
#endif

    if (pivrow [IN] != EMPTY)
    {

    /* the row merge candidate row is not pivrow [IN] */
    freebie [IN] = (pivrow [IN] == prior_pivrow [IN]) && (cdeg1  > 0) ;
    ASSERT (cdeg1  >= 0) ;

    if (!freebie [IN])
    {
        /* include current front in the degree of this row */

        Fcpos = Work->Fcpos ;
        fncols = Work->fncols ;

        Wrpflag = Work->Wrpflag ;

        /* -------------------------------------------------------------- */
        /* construct the pattern of the IN row */
        /* -------------------------------------------------------------- */

#ifndef NDEBUG
        /* check Fcols */
        DEBUG5 (("ROW ASSEMBLE: rdeg "ID"\nREDUCE ROW "ID"\n",
        fncols, pivrow [IN])) ;
        for (j = 0 ; j < fncols ; j++)
        {
        col = Work->Fcols [j] ;
        ASSERT (col >= 0 && col < Work->n_col) ;
        ASSERT (Fcpos [col] >= 0) ;
        }
        if (UMF_debug > 0 || Work->n_col < 1000)
        {
        Int cnt = fncols ;
        for (col = 0 ; col < Work->n_col ; col++)
        {
            if (Fcpos [col] < 0) cnt++ ;
        }
        ASSERT (cnt == Work->n_col) ;
        }
#endif

        rdeg [IN] = fncols ;

        ASSERT (pivrow [IN] >= 0 && pivrow [IN] < Work->n_row) ;
        ASSERT (NON_PIVOTAL_ROW (pivrow [IN])) ;

        /* add the pivot column itself */
        ASSERT (Wrp [pivcol] != Wrpflag) ;
        if (Fcpos [pivcol] < 0)
        {
        DEBUG3 (("Adding pivot col to pivrow [IN] pattern\n")) ;
        if (rdeg [IN] >= max_rdeg)
        {
            /* :: pattern change (in) :: */
            return (UMFPACK_ERROR_different_pattern) ;
        }
        Wrp [pivcol] = Wrpflag ;
        W_i [rdeg [IN]++] = pivcol ;
        }

        tpi = Row_tuples [pivrow [IN]] ;
        if (tpi)
        {
        tp = (Tuple *) (Memory + tpi) ;
        tp1 = tp ;
        tp2 = tp ;
        tpend = tp + Row_tlen [pivrow [IN]] ;
        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 ;
            ncols = ep->ncols ;
            Rows = Cols + ncols ;
            if (Rows [f] == EMPTY)
            {
            continue ;  /* row already assembled */
            }
            ASSERT (pivrow [IN] == Rows [f]) ;

            for (j = 0 ; j < ncols ; j++)
            {
            col = Cols [j] ;
            ASSERT (col >= EMPTY && col < Work->n_col) ;
            if ((col >= 0) && (Wrp [col] != Wrpflag)
                && Fcpos [col] <0)
            {
                ASSERT (NON_PIVOTAL_COL (col)) ;
                if (rdeg [IN] >= max_rdeg)
                {
                /* :: pattern change (rdeg in failure) :: */
                DEBUGm4 (("rdeg [IN] >= max_rdeg failure\n")) ;
                return (UMFPACK_ERROR_different_pattern) ;
                }
                Wrp [col] = Wrpflag ;
                W_i [rdeg [IN]++] = col ;
            }
            }

            *tp2++ = *tp ;  /* leave the tuple in the list */
        }
        Row_tlen [pivrow [IN]] = tp2 - tp1 ;
        }

#ifndef NDEBUG
        DEBUG4 (("Reduced IN row:\n")) ;
        for (j = 0 ; j < fncols ; j++)
        {
        DEBUG6 ((" "ID" "ID" "ID"\n",
            j, Work->Fcols [j], Fcpos [Work->Fcols [j]])) ;
        ASSERT (Fcpos [Work->Fcols [j]] >= 0) ;
        }
        for (j = fncols ; j < rdeg [IN] ; j++)
        {
        DEBUG6 ((" "ID" "ID" "ID"\n", j, W_i [j], Wrp [W_i [j]]));
        ASSERT (W_i [j] >= 0 && W_i [j] < Work->n_col) ;
        ASSERT (Wrp [W_i [j]] == Wrpflag) ;
        }
        /* mark the end of the pattern in case we scan it by mistake */
        /* Note that this means W_i must be of size >= fncols_max + 1 */
        W_i [rdeg [IN]] = EMPTY ;
#endif

        /* rdeg [IN] is now the exact degree of the IN row */

        /* clear Work->Wrp. */
        Work->Wrpflag++ ;
        /* All Wrp [0..n_col] is now < Wrpflag */
    }
    }

#ifndef NDEBUG
    /* check Wrp */
    if (UMF_debug > 0 || Work->n_col < 1000)
    {
    for (i = 0 ; i < Work->n_col ; i++)
    {
        ASSERT (Wrp [i] < Work->Wrpflag) ;
    }
    }
#endif

    /* ---------------------------------------------------------------------- */
    /* construct the candidate row not in the front, if any */
    /* ---------------------------------------------------------------------- */

#ifndef NDEBUG
    DEBUG4 (("pivrow [OUT]: "ID"\n", pivrow [OUT])) ;
    UMF_dump_rowcol (0, Numeric, Work, pivrow [OUT], TRUE) ;
#endif

    /* If this is a candidate row from the row merge set, force it to be */
    /* scanned (ignore prior_pivrow [OUT]). */

    if (pivrow [OUT] != EMPTY)
    {
    freebie [OUT] = (pivrow [OUT] == prior_pivrow [OUT]) && cdeg1  > 0 ;
    ASSERT (cdeg1  >= 0) ;

    if (!freebie [OUT])
    {

        Wrpflag = Work->Wrpflag ;

        /* -------------------------------------------------------------- */
        /* construct the pattern of the row */
        /* -------------------------------------------------------------- */

        rdeg [OUT] = 0 ;

        ASSERT (pivrow [OUT] >= 0 && pivrow [OUT] < Work->n_row) ;
        ASSERT (NON_PIVOTAL_ROW (pivrow [OUT])) ;

        /* add the pivot column itself */
        ASSERT (Wrp [pivcol] != Wrpflag) ;
        DEBUG3 (("Adding pivot col to pivrow [OUT] pattern\n")) ;
        if (rdeg [OUT] >= max_rdeg)
        {
        /* :: pattern change (out) :: */
        return (UMFPACK_ERROR_different_pattern) ;
        }
        Wrp [pivcol] = Wrpflag ;
        W_o [rdeg [OUT]++] = pivcol ;

        tpi = Row_tuples [pivrow [OUT]] ;
        if (tpi)
        {
        tp = (Tuple *) (Memory + tpi) ;
        tp1 = tp ;
        tp2 = tp ;
        tpend = tp + Row_tlen [pivrow [OUT]] ;
        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 ;
            ncols = ep->ncols ;
            Rows = Cols + ncols ;
            if (Rows [f] == EMPTY)
            {
            continue ;  /* row already assembled */
            }
            ASSERT (pivrow [OUT] == Rows [f]) ;

            for (j = 0 ; j < ncols ; j++)
            {
            col = Cols [j] ;
            ASSERT (col >= EMPTY && col < Work->n_col) ;
            if ((col >= 0) && (Wrp [col] != Wrpflag))
            {
                ASSERT (NON_PIVOTAL_COL (col)) ;
                if (rdeg [OUT] >= max_rdeg)
                {
                /* :: pattern change (rdeg out failure) :: */
                DEBUGm4 (("rdeg [OUT] failure\n")) ;
                return (UMFPACK_ERROR_different_pattern) ;
                }
                Wrp [col] = Wrpflag ;
                W_o [rdeg [OUT]++] = col ;
            }
            }
            *tp2++ = *tp ;  /* leave the tuple in the list */
        }
        Row_tlen [pivrow [OUT]] = tp2 - tp1 ;
        }

#ifndef NDEBUG
        DEBUG4 (("Reduced row OUT:\n")) ;
        for (j = 0 ; j < rdeg [OUT] ; j++)
        {
        DEBUG6 ((" "ID" "ID" "ID"\n", j, W_o [j], Wrp [W_o [j]])) ;
        ASSERT (W_o [j] >= 0 && W_o [j] < Work->n_col) ;
        ASSERT (Wrp [W_o [j]] == Wrpflag) ;
        }
        /* mark the end of the pattern in case we scan it by mistake */
        /* Note that this means W_o must be of size >= fncols_max + 1 */
        W_o [rdeg [OUT]] = EMPTY ;
#endif

        /* rdeg [OUT] is now the exact degree of the row */

        /* clear Work->Wrp. */
        Work->Wrpflag++ ;
        /* All Wrp [0..n] is now < Wrpflag */

    }

    }
    DEBUG5 (("Row_search end ] \n")) ;

#ifndef NDEBUG
    /* check Wrp */
    if (UMF_debug > 0 || Work->n_col < 1000)
    {
    for (i = 0 ; i < Work->n_col ; i++)
    {
        ASSERT (Wrp [i] < Work->Wrpflag) ;
    }
    }
#endif

    return (UMFPACK_OK) ;
}