The R Project SVN R-packages

/* ========================================================================== */
/* === UMF_dump ============================================================= */
/* ========================================================================== */

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

/* These routines, and external variables, are used only when debugging. */
/* If debugging is disabled (for normal operation) then this entire file */
/* becomes empty */

#include "umf_internal.h"

#ifndef NDEBUG

/* These global debugging variables and arrays do not exist if debugging */
/* is disabled at compile time (which is the default). */
GLOBAL Int UMF_debug = -999 ;
GLOBAL Int UMF_allocfail = FALSE ;
GLOBAL double UMF_gprob = -1.0 ;

/* static debugging arrays used only in UMF_dump_rowcol */
PRIVATE Int UMF_DBflag = 0 ;
PRIVATE Int UMF_DBpacked [UMF_DBMAX+1] ;
PRIVATE Int UMF_DBscatter [UMF_DBMAX+1] ;

/* ========================================================================== */
/* === UMF_DBinit =========================================================== */
/* ========================================================================== */

/* clear the debugging arrays */

PRIVATE void UMF_DBinit
(
    void
)
{
    Int i ;

    /* Int_MAX is defined in umfpack.h */
    if (UMF_DBflag < 1 || UMF_DBflag == Int_MAX)
    {
    /* clear the debugging arrays */
    UMF_DBflag = 0 ;
    for (i = 0 ; i <= UMF_DBMAX ; i++)
    {
        UMF_DBscatter [i] = 0 ;
        UMF_DBpacked  [i] = 0 ;
    }
    }

    UMF_DBflag++ ;

    /* UMF_DBflag > UMF_DBscatter [0...UMF_DBmax] is now true */
}

/* ========================================================================== */
/* === UMF_dump_dense ======================================================= */
/* ========================================================================== */

GLOBAL void UMF_dump_dense
(
    Entry *C,
    Int dim,
    Int m,
    Int n
)
{

    /* dump C [1..m,1..n], with column dimenstion dim */
    Int i, j;

    if (UMF_debug < 7) return ;
    if (C == (Entry *) NULL)
    {
    DEBUG7 (("No dense matrix allocated\n")) ;
    return ;
    }
    DEBUG8 ((" dimension= "ID" rows= "ID" cols= "ID"\n", dim, m, n)) ;

    for (i = 0 ; i < m ; i++)
    {
    DEBUG9 ((ID": ", i)) ;
    for (j = 0 ; j < n ; j++)
    {
        EDEBUG9 (C [i+j*dim]) ;
        if (j % 6 == 5) DEBUG9 (("\n     ")) ;
    }
    DEBUG9 (("\n")) ;
    }

    for (i = 0 ; i < m ; i++)
    {
    for (j = 0 ; j < n ; j++)
    {
        if (IS_ZERO (C [i+j*dim]))
        {
        DEBUG8 ((".")) ;
        }
        else
        {
        DEBUG8 (("X")) ;
        }
    }
    DEBUG8 (("\n")) ;
    }
}

/* ========================================================================== */
/* === UMF_dump_element ===================================================== */
/* ========================================================================== */

GLOBAL void UMF_dump_element
(
    NumericType *Numeric,
    WorkType *Work,
    Int e,
    Int clean
)
{

    Int i, j, k, *Rows, *Cols, nrows, ncols, *E, row, col,
    *Row_degree, *Col_degree ;
    Entry *C ;
    Element *ep ;
    Unit *p ;

    if (UMF_debug < 7) return ;

    if (e == 0)
    {
    UMF_dump_current_front (Numeric, Work, FALSE) ;
    return ;
    }

    DEBUG7 (("\n====================ELEMENT: "ID" ", e)) ;
    if (!Numeric || !Work || !Numeric->Memory)
    {
    DEBUG7 ((" No Numeric, Work\n")) ;
    return ;
    }
    DEBUG7 ((" nel: "ID" of "ID, e, Work->nel)) ;
    E = Work->E ;
    if (!E)
    {
    DEBUG7 ((" No elements\n")) ;
    return ;
    }
    if (e < 0 || e > Work->nel)
    {
    DEBUG7 (("e out of range!\n")) ;
    return ;
    }
    if (!E [e])
    {
    DEBUG7 ((" deallocated\n")) ;
    return ;
    }
    DEBUG7 (("\n")) ;
    Col_degree = Numeric->Cperm ;
    Row_degree = Numeric->Rperm ;

    p = Numeric->Memory + E [e] ;
    DEBUG7 (("ep "ID"\n", (Int) (p-Numeric->Memory))) ;
    GET_ELEMENT (ep, p, Cols, Rows, ncols, nrows, C) ;
    DEBUG7 (("nrows "ID" nrowsleft "ID"\n", nrows, ep->nrowsleft)) ;
    DEBUG7 (("ncols "ID" ncolsleft "ID"\n", ncols, ep->ncolsleft)) ;
    DEBUG7 (("cdeg-cdeg0 "ID" rdeg-rdeg0 "ID" next "ID"\n",
    ep->cdeg - Work->cdeg0, ep->rdeg - Work->rdeg0, ep->next)) ;

    DEBUG8 (("rows: ")) ;
    k = 0 ;
    for (i = 0 ; i < ep->nrows ; i++)
    {
    row = Rows [i] ;
    if (row >= 0)
    {
        DEBUG8 ((" "ID, row)) ;
        ASSERT (row < Work->n_row) ;
        if ((k++ % 10) == 9) DEBUG8 (("\n")) ;
        ASSERT (IMPLIES (clean, NON_PIVOTAL_ROW (row))) ;
    }
    }

    DEBUG8 (("\ncols: ")) ;
    k = 0 ;
    for (j = 0 ; j < ep->ncols ; j++)
    {
    col = Cols [j] ;
    if (col >= 0)
    {
        DEBUG8 ((" "ID, col)) ;
        ASSERT (col < Work->n_col) ;
        if ((k++ % 10) == 9) DEBUG8 (("\n")) ;
        ASSERT (IMPLIES (clean, NON_PIVOTAL_COL (col))) ;
    }
    }

    DEBUG8 (("\nvalues:\n")) ;
    if (UMF_debug >= 9)
    {
    for (i = 0 ; i < ep->nrows ; i++)
    {
        row = Rows [i] ;
        if (row >= 0)
        {
        DEBUG9 ((ID": ", row)) ;
        k = 0 ;
        for (j = 0 ; j < ep->ncols ; j++)
        {
            col = Cols [j] ;
            if (col >= 0)
            {
            EDEBUG9 (C [i+j*ep->nrows]) ;
            if (k++ % 6 == 5) DEBUG9 (("\n     ")) ;
            }
        }
        DEBUG9 (("\n")) ;
        }
    }
    }

    DEBUG7 (("====================\n")) ;
}


/* ========================================================================== */
/* === UMF_dump_rowcol ====================================================== */
/* ========================================================================== */

/* dump a row or a column, from one or more memory spaces */
/* return exact degree */

GLOBAL void UMF_dump_rowcol
(
    Int dumpwhich,      /* 0 for row, 1 for column */
    NumericType *Numeric,
    WorkType *Work,
    Int dumpindex,      /* row or column index to dump */
    Int check_degree    /* true if degree is to be checked */
)
{
    Entry value ;
    Entry *C ;
    Int f, nrows, j, jj, len, e, deg, index, n_row, n_col, *Cols, *Rows, nn,
    dumpdeg, ncols, preve, *E, tpi, *Pattern, approx_deg, not_in_use ;
    Tuple *tp, *tend ;
    Element *ep ;
    Int *Row_tuples, *Row_degree, *Row_tlen ;
    Int *Col_tuples, *Col_degree, *Col_tlen ;
    Unit *p ;
    Int is_there ;

    /* clear the debugging arrays */
    UMF_DBinit () ;

    if (dumpwhich == 0)
    {
    DEBUG7 (("\n====================ROW: "ID, dumpindex)) ;
    }
    else
    {
    DEBUG7 (("\n====================COL: "ID, dumpindex)) ;
    }

    if (dumpindex == EMPTY)
    {
    DEBUG7 ((" (EMPTY)\n")) ;
    return ;
    }

    deg = 0 ;
    approx_deg = 0 ;

    if (!Numeric || !Work)
    {
    DEBUG7 ((" No Numeric, Work\n")) ;
    return ;
    }

    n_row = Work->n_row ;
    n_col = Work->n_col ;
    nn = MAX (n_row, n_col) ;
    E = Work->E ;

    Col_degree = Numeric->Cperm ;
    Row_degree = Numeric->Rperm ;

    Row_tuples = Numeric->Uip ;
    Row_tlen   = Numeric->Uilen ;
    Col_tuples = Numeric->Lip ;
    Col_tlen   = Numeric->Lilen ;

    if (!E
    || !Row_tuples || !Row_degree || !Row_tlen
    || !Col_tuples || !Col_degree || !Col_tlen)
    {
        DEBUG7 ((" No E, Rows, Cols\n")) ;
        return ;
    }

    if (dumpwhich == 0)
    {
        /* dump a row */
        ASSERT (dumpindex >= 0 && dumpindex < n_row) ;
        if (!NON_PIVOTAL_ROW (dumpindex))
        {
        DEBUG7 ((" Pivotal\n")) ;
        return ;
        }
        len = Row_tlen [dumpindex] ;
        dumpdeg = Row_degree [dumpindex] ;
        tpi = Row_tuples [dumpindex] ;
    }
    else
    {
        /* dump a column */
        ASSERT (dumpindex >= 0 && dumpindex < n_col) ;
        if (!NON_PIVOTAL_COL (dumpindex))
        {
        DEBUG7 ((" Pivotal\n")) ;
        return ;
        }
        len = Col_tlen [dumpindex] ;
        dumpdeg = Col_degree [dumpindex] ;
        tpi = Col_tuples [dumpindex] ;
    }

    p = Numeric->Memory + tpi ;
    tp = (Tuple *) p ;
    if (!tpi)
    {
        DEBUG7 ((" Nonpivotal, No tuple list tuples "ID" tlen "ID"\n",
        tpi, len)) ;
        return ;
    }
    ASSERT (p >= Numeric->Memory + Numeric->itail) ;
    ASSERT (p <  Numeric->Memory + Numeric->size) ;

    DEBUG7 ((" degree: "ID" len: "ID"\n", dumpdeg, len)) ;
    not_in_use = (p-1)->header.size - UNITS (Tuple, len) ;
    DEBUG7 ((" Tuple list: p+1: "ID" size: "ID" units, "ID" not in use\n",
        (Int) (p-Numeric->Memory), (p-1)->header.size, not_in_use)) ;
    ASSERT (not_in_use >= 0) ;
    tend = tp + len ;
    preve = 0 ;
    for ( ; tp < tend ; tp++)
    {
        /* row/col of element e, offset is f: */
        /* DEBUG8 (("    (tp="ID")\n", tp)) ; */
        e = tp->e ;
        f = tp->f ;
        DEBUG8 (("    (e="ID", f="ID")\n", e, f)) ;
        ASSERT (e > 0 && e <= Work->nel) ;
        /* dump the pattern and values */
        if (E [e])
        {
        p = Numeric->Memory + E [e] ;
        GET_ELEMENT (ep, p, Cols, Rows, ncols, nrows, C) ;
        if (dumpwhich == 0)
        {
            Pattern = Cols ;
            jj = ep->ncols ;
            is_there = Rows [f] >= 0 ;
            if (is_there) approx_deg += ep->ncolsleft ;
        }
        else
        {
            Pattern = Rows ;
            jj = ep->nrows ;
            is_there = Cols [f] >= 0 ;
            if (is_there) approx_deg += ep->nrowsleft ;
        }
        if (!is_there)
        {
            DEBUG8 (("\t\tnot present\n")) ;
        }
        else
        {
            for (j = 0 ; j < jj ; j++)
            {
            index = Pattern [j] ;
            value =
                C [ (dumpwhich == 0) ? (f+nrows*j) : (j+nrows*f) ] ;
            if (index >= 0)
            {
                DEBUG8 (("\t\t"ID":", index)) ;
                EDEBUG8 (value) ;
                DEBUG8 (("\n")) ;
                if (dumpwhich == 0)
                {
                /* col must be in the range 0..n_col-1 */
                ASSERT (index < n_col) ;
                }
                else
                {
                /* row must be in the range 0..n_row-1 */
                ASSERT (index < n_row) ;
                }

                if (nn <= UMF_DBMAX)
                {
                if (UMF_DBscatter [index] != UMF_DBflag)
                {
                    UMF_DBpacked [deg++] = index ;
                    UMF_DBscatter [index] = UMF_DBflag ;
                }
                }
            }
            }
        }
        /* the (e,f) tuples should be in order of their creation */
        /* this means that garbage collection will not jumble them */
        ASSERT (preve < e) ;
        preve = e ;
        }
        else
        {
        DEBUG8 (("\t\tdeallocated\n")) ;
        }
    }

    if (nn <= UMF_DBMAX)
    {
    if (deg > 0)
    {
        DEBUG7 ((" Assembled, actual deg: "ID" : ", deg)) ;
        for (j = 0 ; j < deg ; j++)
        {
        index = UMF_DBpacked [j] ;
        DEBUG8 ((ID" ", index)) ;
        if (j % 20 == 19) DEBUG8 (("\n ")) ;
        ASSERT (UMF_DBscatter [index] == UMF_DBflag) ;
        }
        DEBUG7 (("\n")) ;
    }
    }

    /* Col_degree is not maintained when fixQ is true */
    if (check_degree)
    {
    DEBUG8 (("  approx_deg "ID"  dumpdeg "ID"\n", approx_deg, dumpdeg)) ;
    ASSERT (approx_deg == dumpdeg) ;
    }

    DEBUG7 (("====================\n")) ;

    /* deg is now the exact degree */
    /* if nn <= UMF_DBMAX, then UMF_DBscatter [i] == UMF_DBflag for every i */
    /* in the row or col, and != UMF_DBflag if not */

    return ;
}


/* ========================================================================== */
/* === UMF_dump_matrix ====================================================== */
/* ========================================================================== */

GLOBAL void UMF_dump_matrix
(
    NumericType *Numeric,
    WorkType *Work,
    Int check_degree
)
{

    Int e, row, col, intfrag, frag, n_row, n_col, *E, fullsize, actualsize ;
    Element *ep ;
    Unit *p ;

    DEBUG6 (("=================================================== MATRIX:\n")) ;
    if (!Numeric || !Work)
    {
    DEBUG6 (("No Numeric or Work allocated\n")) ;
    return ;
    }
    if (!Numeric->Memory)
    {
    DEBUG6 (("No Numeric->Memory\n")) ;
    return ;
    }

    n_row = Work->n_row ;
    n_col = Work->n_col ;
    DEBUG6 (("n_row "ID" n_col "ID" nz "ID"\n", n_row, n_col, Work->nz)) ;
    DEBUG6 (("============================ ELEMENTS: "ID" \n", Work->nel)) ;
    intfrag = 0 ;
    E = Work->E ;
    if (!E)
    {
        DEBUG6 (("No elements allocated\n")) ;
    }
    else
    {
        for (e = 0 ; e <= Work->nel ; e++)
        {
        UMF_dump_element (Numeric, Work, e, FALSE) ;
        if (e > 0 && E [e])
        {
            p = Numeric->Memory + E [e] ;
            ep = (Element *) p ;
            ASSERT (ep->nrowsleft > 0 || ep->ncolsleft > 0) ;
            fullsize = GET_BLOCK_SIZE (p) ;
            actualsize = GET_ELEMENT_SIZE (ep->nrowsleft,ep->ncolsleft);
            frag =  fullsize - actualsize ;
            intfrag += frag ;
            DEBUG7 (("dump el: "ID", full "ID" actual "ID" frag: "ID
            " intfrag: "ID"\n", e, fullsize, actualsize, frag,
            intfrag)) ;
        }
        }
    }

    DEBUG6 (("CURRENT INTERNAL FRAG in elements: "ID" \n", intfrag)) ;



    DEBUG6 (("======================================== ROWS: "ID"\n", n_row)) ;
    UMF_debug -= 2 ;
    for (row = 0 ; row < n_row ; row++)
    {
    UMF_dump_rowcol (0, Numeric, Work, row, check_degree) ;
    }
    UMF_debug += 2 ;
    DEBUG6 (("======================================== COLS: "ID"\n", n_col)) ;
    UMF_debug -= 2 ;
    for (col = 0 ; col < n_col ; col++)
    {
    UMF_dump_rowcol (1, Numeric, Work, col, FALSE) ;
    }
    UMF_debug += 2 ;
    DEBUG6 (("============================================= END OF MATRIX:\n"));
}


/* ========================================================================== */
/* === UMF_dump_current_front =============================================== */
/* ========================================================================== */

GLOBAL void UMF_dump_current_front
(
    NumericType *Numeric,
    WorkType *Work,
    Int check
)
{

    Entry *Flublock, *Flblock, *Fublock, *Fcblock ;
    Int fnrows_max, fncols_max, fnrows, fncols, fnpiv, *Frows, *Fcols,
    i, j, *Fcpos, *Frpos, fnr_curr, fnc_curr, *E ;
    if (!Work) return ;
    DEBUG7 (("\n\n========CURRENT FRONTAL MATRIX:\n")) ;

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

    Frows = Work->Frows ;
    Fcols = Work->Fcols ;
    Frpos = Work->Frpos ;
    Fcpos = Work->Fcpos ;
    fnrows_max = Work->fnrows_max ;
    fncols_max = Work->fncols_max ;
    fnr_curr = Work->fnr_curr ;
    fnc_curr = Work->fnc_curr ;
    fnrows = Work->fnrows ;
    fncols = Work->fncols ;
    fnpiv = Work->fnpiv ;
    E = Work->E ;

    DEBUG6 (("=== fnpiv= "ID"\n", fnpiv)) ;
    DEBUG6 (("fnrows_max      fncols_max "ID" "ID"\n",fnrows_max, fncols_max)) ;
    DEBUG6 (("fnr_curr        fnc_curr   "ID" "ID"\n",fnr_curr,   fnc_curr)) ;
    DEBUG6 (("fnrows          fncols     "ID" "ID"\n",fnrows,     fncols)) ;
    ASSERT ((fnr_curr % 2 == 1) || fnr_curr == 0) ;
    DEBUG6 (("Pivot row pattern:\n")) ;
    for (j = 0 ; j < fncols ; j++)
    {
    DEBUG7 ((ID" "ID" "ID" %d\n", j, Fcols [j], Fcpos [Fcols [j]],
        j < fncols)) ;
    if (check)
    {
        ASSERT (Fcols [j] >= 0 && Fcols [j] < Work->n_col) ;
        ASSERT (Fcpos [Fcols [j]] == j * fnr_curr) ;
    }
    }
    DEBUG6 (("Pivot col pattern:\n")) ;
    for (i = 0 ; i < fnrows ; i++)
    {
    DEBUG7 ((ID" "ID" "ID" %d\n", i, Frows [i], Frpos [Frows [i]],
        i < fnrows)) ;
    if (check)
    {
        ASSERT (Frows [i] >= 0 && Frows [i] < Work->n_row) ;
        ASSERT (Frpos [Frows [i]] == i) ;
    }
    }
    if (UMF_debug < 7) return ;

    if (!E [0])
    {
    DEBUG6 (("current front not allocated\n")) ;
    ASSERT (!Work->Flublock) ;
    return ;
    }

    ASSERT (Work->Flublock == (Entry *) (Numeric->Memory + E [0])) ;
    DEBUG7 (("C  block: ")) ;
    UMF_dump_dense (Fcblock,  fnr_curr, fnrows, fncols) ;
    DEBUG7 (("L  block: ")) ;
    UMF_dump_dense (Flblock,  fnr_curr, fnrows, fnpiv) ;
    DEBUG7 (("U' block: ")) ;
    UMF_dump_dense (Fublock,  fnc_curr, fncols, fnpiv) ;
    DEBUG7 (("LU block: ")) ;
    UMF_dump_dense (Flublock, Work->nb, fnpiv, fnpiv) ;
    if (fnpiv > 0)
    {
    DEBUG7 (("Pivot entry: ")) ;
    EDEBUG7 (Flublock [(fnpiv-1)+(fnpiv-1)*Work->nb]) ;
    DEBUG7 (("\n")) ;
    }
}

/* ========================================================================== */
/* === UMF_dump_lu ========================================================== */
/* ========================================================================== */

GLOBAL void UMF_dump_lu
(
    NumericType *Numeric
)
{
    Int i, n_row, n_col, *Cperm, *Rperm ;

    DEBUG6 (("=============================================== LU factors:\n")) ;
    if (!Numeric)
    {
    DEBUG6 (("No LU factors allocated\n")) ;
    return ;
    }
    n_row = Numeric->n_row ;
    n_col = Numeric->n_col ;
    DEBUG6 (("n_row: "ID" n_col: "ID"\n", n_row, n_col)) ;
    DEBUG6 (("nLentries: "ID" nUentries: "ID"\n",
    Numeric->nLentries, Numeric->nUentries)) ;

    if (Numeric->Cperm)
    {
    Cperm = Numeric->Cperm ;
    DEBUG7 (("Column permutations: (new: old)\n")) ;
    for (i = 0 ; i < n_col ; i++)
    {
        if (Cperm [i] != EMPTY)
        {
        DEBUG7 ((ID": "ID"\n", i, Cperm [i])) ;
        }
    }
    }
    else
    {
    DEBUG7 (("No Numeric->Cperm allocatated\n")) ;
    }

    if (Numeric->Rperm)
    {
    Rperm = Numeric->Rperm ;
    DEBUG7 (("row permutations: (new: old)\n")) ;
    for (i = 0 ; i < n_row ; i++)
    {
        if (Rperm [i] != EMPTY)
        {
        DEBUG7 ((ID": "ID"\n", i, Rperm [i])) ;
        }
    }
    }
    else
    {
    DEBUG7 (("No Numeric->Rperm allocatated\n")) ;
    }

    DEBUG6 (("========================================= END OF LU factors:\n"));
}


/* ========================================================================== */
/* === UMF_dump_memory ====================================================== */
/* ========================================================================== */

GLOBAL void UMF_dump_memory
(
    NumericType *Numeric
)
{

    Unit *p ;
    Int prevsize, s ;
    Int found ;

    if (!Numeric)
    {
    DEBUG6 (("No memory space S allocated\n")) ;
    return ;
    }

    DEBUG6 (("\n ============================================== MEMORY:\n")) ;
    if (!Numeric || !Numeric->Memory)
    {
    DEBUG6 (("No memory space Numeric allocated\n")) ;
    return ;
    }

    DEBUG6 (("S: "ID"\n", (Int) Numeric)) ;
    DEBUG6 (("S->ihead           : "ID"\n", Numeric->ihead)) ;
    DEBUG6 (("S->itail           : "ID"\n", Numeric->itail)) ;
    DEBUG6 (("S->size            : "ID"\n", Numeric->size)) ;
    DEBUG6 (("S->ngarbage        : "ID"\n", Numeric->ngarbage)) ;
    DEBUG6 (("S->nrealloc        : "ID"\n", Numeric->nrealloc)) ;
    DEBUG6 (("   in use at head           : "ID"\n", Numeric->ihead)) ;
    DEBUG6 (("   free space               : "ID"\n",
    Numeric->itail - Numeric->ihead)) ;
    DEBUG6 (("   blocks in use at tail    : "ID"\n",
    Numeric->size - Numeric->itail)) ;
    DEBUG6 (("   total in use             : "ID"\n",
    Numeric->size - (Numeric->itail - Numeric->ihead))) ;

    prevsize = 0 ;
    found = FALSE ;

    ASSERT (0 <= Numeric->ihead) ;
    ASSERT (Numeric->ihead <= Numeric->itail) ;
    ASSERT (Numeric->itail <= Numeric->size) ;

    p = Numeric->Memory + Numeric->itail ;

    while (p < Numeric->Memory + Numeric->size)
    {
    DEBUG8 (("p: "ID" p+1: "ID" prevsize: "ID" size: "ID,
        (Int) (p-Numeric->Memory), (Int) (p+1-Numeric->Memory),
        p->header.prevsize, p->header.size)) ;
    if (p->header.size < 0)
    {
        DEBUG8 ((" free")) ;
    }

    if (p == Numeric->Memory + Numeric->itail)
    {
        ASSERT (p->header.prevsize == 0) ;
    }
    else
    {
        ASSERT (p->header.prevsize > 0) ;
    }

    ASSERT (p->header.size != 0) ;
    s = prevsize >= 0 ? prevsize : -prevsize ;
    ASSERT (p->header.prevsize == s) ;
    /* no adjacent free blocks */
    ASSERT (p->header.size > 0 || prevsize > 0) ;
    if (Numeric->ibig != EMPTY)
    {
        if (p == Numeric->Memory + Numeric->ibig)
        {
        ASSERT (p->header.size < 0) ;
        DEBUG8 ((" <===== Numeric->ibig")) ;
        found = TRUE ;
        }
    }
    s = p->header.size ;
    prevsize = s ;
    s = s >= 0 ? s : -s ;
    p = p + 1 + s ;
    DEBUG8 (("\n")) ;
    }

    ASSERT (p == Numeric->Memory + Numeric->size) ;
    ASSERT (IMPLIES (Numeric->ibig != EMPTY, found)) ;
    DEBUG6 (("============================================= END OF MEMORY:\n"));

}


/* ========================================================================== */
/* === UMF_dump_packed_memory =============================================== */
/* ========================================================================== */

GLOBAL void UMF_dump_packed_memory
(
    NumericType *Numeric,
    WorkType *Work
)
{
    Unit *p, *p3 ;
    Int prevsize, col, row, *Rows, *Cols, ncols, nrows, k, esize,
    *Row_tuples, *Row_degree, *Col_tuples, *Col_degree ;
    Entry *C ;
    Element *ep ;

    Col_degree = Numeric->Cperm ;   /* for NON_PIVOTAL_COL macro */
    Row_degree = Numeric->Rperm ;   /* for NON_PIVOTAL_ROW macro */
    Row_tuples = Numeric->Uip ;
    Col_tuples = Numeric->Lip ;

    DEBUG6 (("============================================ PACKED MEMORY:\n")) ;
    if (!Numeric || !Numeric->Memory)
    {
    DEBUG6 (("No memory space S allocated\n")) ;
    return ;
    }
    DEBUG6 (("S: "ID"\n", (Int) Numeric)) ;
    DEBUG6 (("S->ihead           : "ID"\n", Numeric->ihead)) ;
    DEBUG6 (("S->itail           : "ID"\n", Numeric->itail)) ;
    DEBUG6 (("S->size            : "ID"\n", Numeric->size)) ;
    DEBUG6 (("S->ngarbage        : "ID"\n", Numeric->ngarbage)) ;
    DEBUG6 (("S->nrealloc        : "ID"\n", Numeric->nrealloc)) ;
    DEBUG6 (("   in use at head           : "ID"\n", Numeric->ihead)) ;
    DEBUG6 (("   free space               : "ID"\n",
    Numeric->itail - Numeric->ihead)) ;
    DEBUG6 (("   blocks in use at tail    : "ID"\n",
    Numeric->size - Numeric->itail)) ;
    DEBUG6 (("   total in use             : "ID"\n",
    Numeric->size - (Numeric->itail - Numeric->ihead))) ;

    ASSERT (0 <= Numeric->ihead) ;
    ASSERT (Numeric->ihead <= Numeric->itail) ;
    ASSERT (Numeric->itail <= Numeric->size) ;

    for (row = 0 ; row < Work->n_row ; row++)
    {
    ASSERT (IMPLIES (NON_PIVOTAL_ROW (row), !Row_tuples [row])) ;
    }
    for (col = 0 ; col < Work->n_col ; col++)
    {
    ASSERT (IMPLIES (NON_PIVOTAL_COL (col), !Col_tuples [col])) ;
    }

    prevsize = 0 ;
    p = Numeric->Memory + Numeric->itail ;
    while (p < Numeric->Memory + Numeric->size)
    {
    DEBUG9 (("====================\n")) ;
    DEBUG7 (("p: "ID" p+1: "ID" prevsize: "ID" size: "ID"\n",
        (Int) (p-Numeric->Memory), (Int) (p+1-Numeric->Memory),
        p->header.prevsize, p->header.size)) ;
    ASSERT (p->header.size > 0) ;

    if (p == Numeric->Memory + Numeric->itail)
    {
        ASSERT (p->header.prevsize == 0) ;
    }
    else
    {
        ASSERT (p->header.prevsize > 0) ;
    }

    ASSERT (p->header.prevsize == prevsize) ;
    prevsize = p->header.size ;

    if (p != Numeric->Memory + Numeric->size - 2)
    {

        p3 = p + 1 ;
        if (p3 == Numeric->Memory + Work->E [0])
        {
        /* this is the current frontal matrix */
        UMF_dump_current_front (Numeric, Work, FALSE) ;
        }
        else
        {

        /* this is a packed element */
        GET_ELEMENT (ep, p3, Cols, Rows, ncols, nrows, C) ;
        DEBUG9 (("ep "ID"\n nrows "ID" ncols "ID"\n",
            (Int) ((p+1)-Numeric->Memory), ep->nrows, ep->ncols)) ;
        DEBUG9 (("rows:")) ;
        for (k = 0 ; k < ep->nrows; k++)
        {
            row = Rows [k] ;
            DEBUG9 ((" "ID, row)) ;
            ASSERT (row >= 0 && row <= Work->n_row) ;
            if ((k % 10) == 9) DEBUG9 (("\n")) ;
        }
        DEBUG9 (("\ncols:")) ;
        for (k = 0 ; k < ep->ncols; k++)
        {
            col = Cols [k] ;
            DEBUG9 ((" "ID, col)) ;
            ASSERT (col >= 0 && col <= Work->n_col) ;
            if ((k % 10) == 9) DEBUG9 (("\n")) ;
        }
        DEBUG9 (("\nvalues: ")) ;
        if (UMF_debug >= 9)
        {
            UMF_dump_dense (C, ep->nrows, ep->nrows, ep->ncols) ;
        }
        esize = GET_ELEMENT_SIZE (ep->nrows, ep->ncols) ;
        DEBUG9 (("esize: "ID"\n", esize)) ;
        ASSERT (esize <= p->header.size) ;
        }

    }
    else
    {
        /* this is the final marker block */
        ASSERT (p->header.size == 1) ;
    }
    p = p + 1 + p->header.size ;
    }

    ASSERT (Numeric->ibig == EMPTY) ;
    ASSERT (p == Numeric->Memory + Numeric->size) ;
    DEBUG6 (("======================================END OF PACKED MEMORY:\n")) ;

}

/* ========================================================================== */
/* === UMF_dump_col_matrix ================================================== */
/* ========================================================================== */

/* This code is the same for real or complex matrices. */

GLOBAL void UMF_dump_col_matrix
(
    const double Ax [ ],    /* Ax [0..nz-1]: real values, in column order */
#ifdef COMPLEX
    const double Az [ ],    /* Az [0..nz-1]: imag values, in column order */
#endif
    const Int Ai [ ],       /* Ai [0..nz-1]: row indices, in column order */
    const Int Ap [ ],       /* Ap [0..n_col]: column pointers */
    Int n_row,          /* number of rows of A */
    Int n_col,          /* number of columns of A */
    Int nz          /* number of entries */
)
{
    Int col, p, p1, p2, row ;
#ifdef COMPLEX
    Int split = SPLIT (Az) ;
#endif

    if (!Ai || !Ap) return ;
    DEBUG6 (("============================================ COLUMN FORM:\n")) ;

    ASSERT (n_col >= 0) ;
    nz = Ap [n_col] ;
    DEBUG2 (("UMF_dump_col:  nz "ID"\n", nz)) ;
    DEBUG2 (("n_row "ID"  \n", n_row)) ;
    DEBUG2 (("n_col "ID"  \n", n_col)) ;

    DEBUG6 ((" n_row = "ID", n_col ="ID" nz = "ID" Ap [0] "ID", Ap [n] "ID"\n",
    n_row, n_col, nz, Ap [0], Ap [n_col])) ;
    ASSERT (Ap [0] == 0) ;
    ASSERT (Ap [n_col] == nz) ;
    for (col = 0 ; col < n_col ; col++)
    {
    p1 = Ap [col] ;
    p2 = Ap [col+1] ;
    DEBUG6 (("col: "ID", length "ID"\n", col, p2 - p1)) ;
    ASSERT (p2 >= p1) ;
    for (p = p1 ; p < p2 ; p++)
    {
        row = Ai [p] ;
        ASSERT (row >= 0 && row < n_row) ;
        DEBUG6 (("\t"ID" ", row)) ;
        if (Ax != (double *) NULL)
        {
#ifdef COMPLEX
        if (split)
        {
            DEBUG6 ((" (%e+%ei) ", Ax [p], Az [p])) ;
        }
        else
        {
            DEBUG6 ((" (%e+%ei) ", Ax [2*p], Ax [2*p+1])) ;
        }
#else
        DEBUG6 ((" %e", Ax [p])) ;
#endif
        }
        DEBUG6 (("\n")) ;
    }
    }
    DEBUG6 (("========================================== COLUMN FORM done\n")) ;
}


/* ========================================================================== */
/* === UMF_dump_chain ======================================================= */
/* ========================================================================== */

GLOBAL void UMF_dump_chain
(
    Int frontid,
    Int Front_parent [ ],
    Int Front_npivcol [ ],
    Int Front_nrows [ ],
    Int Front_ncols [ ],
    Int nfr
)
{
    Int i, len = 0 ;

    /* print a list of contiguous parents */
    i = frontid ;
    ASSERT (Front_parent [i] == EMPTY ||
    (Front_parent [i] > i && Front_parent [i] < nfr)) ;

    len++ ;
    DEBUG3 (("Chain:\n  "ID" ["ID","ID"]("ID"-by-"ID")\n", i,
        Front_npivcol [i],
        MIN (Front_npivcol [i], Front_nrows [i]),
        Front_nrows [i],
        Front_ncols [i])) ;

    for (i = frontid ; i < nfr ; i++)
    {
    ASSERT (Front_parent [i] == EMPTY ||
    (Front_parent [i] > i && Front_parent [i] < nfr)) ;
    if (Front_parent [i] == i+1)
    {
        len++ ;
        DEBUG3 (("\t"ID" ["ID","ID"]("ID"-by-"ID")\n", i+1,
        Front_npivcol [i+1],
        MIN (Front_npivcol [i+1], Front_nrows [i+1]),
        Front_nrows [i+1],
        Front_ncols [i+1])) ;
    }
    else
    {
        DEBUG2 (("Length of chain: "ID"\n", len)) ;
        return ;
    }
    }
}


/* ========================================================================== */
/* === UMF_dump_start ======================================================= */
/* ========================================================================== */

GLOBAL void UMF_dump_start
(
    void
)
{
    FILE *ff ;

    AMD_debug_init ("from umfpack") ;

    /* get the debug print level from the "debug.umf" file, if it exists */
    UMF_debug = -999 ;
    ff = fopen ("debug.umf", "r") ;
    if (ff)
    {
    (void) fscanf (ff, ID, &UMF_debug) ;
    (void) fclose (ff) ;
    }

    DEBUG0 (("umfpack: debug version (SLOW!) ")) ;

    DEBUG0 ((" BLAS: ")) ;

#if defined (USE_NO_BLAS)
    DEBUG0 (("none.")) ;
#elif defined (USE_C_BLAS)
    DEBUG0 (("C-BLAS.")) ;
#elif defined (USE_MATLAB_BLAS)
    DEBUG0 (("built-in MATLAB BLAS.")) ;
#elif defined (USE_SUNPERF_BLAS)
    DEBUG0 (("Sun Performance Library BLAS.")) ;
#elif defined (USE_SCSL_BLAS)
    DEBUG0 (("SGI SCSL BLAS.")) ;
#elif defined (USE_FORTRAN_BLAS)
    DEBUG0 (("Fortran BLAS.")) ;
#endif

    DEBUG0 ((" MATLAB: ")) ;
#ifdef MATLAB_MEX_FILE
    DEBUG0 (("mexFunction.\n")) ;
#else
#ifdef MATHWORKS
    DEBUG0 (("yes.\n")) ;
#else
    DEBUG0 (("no.\n")) ;
#endif
#endif

    UMF_gprob = -1.0 ;
    ff = fopen ("gprob.umf", "r") ;
    if (ff)
    {
    (void) fscanf (ff, "%lg", &UMF_gprob) ;
    (void) fclose (ff) ;
    srand (1) ; /* restart the random number generator */
    }

    if (UMF_gprob > 1.0) UMF_gprob = 1.0 ;
    DEBUG1 (("factor: UMF_gprob: %e UMF_debug "ID"\n", UMF_gprob, UMF_debug)) ;

    DEBUG2 (("sizeof: (bytes / int / Units) \n")) ;
    DEBUG2 (("sizeof (Int)           %u %u %u\n",
    sizeof (Int), sizeof (Int) / sizeof (int), UNITS (Int, 1) )) ;
    DEBUG2 (("sizeof (int)           %u %u %u\n",
    sizeof (int), sizeof (int) / sizeof (int), UNITS (int, 1) )) ;
    DEBUG2 (("sizeof (size_t)        %u %u %u\n",
    sizeof (size_t), sizeof (size_t) / sizeof (size_t), UNITS (size_t, 1) )) ;
    DEBUG2 (("sizeof (long)          %u %u %u\n",
    sizeof (long), sizeof (long) / sizeof (long), UNITS (long, 1) )) ;
    DEBUG2 (("sizeof (double)        %u %u %u\n",
    sizeof (double), sizeof (double) / sizeof (int), UNITS (double, 1) )) ;
    DEBUG2 (("sizeof (Unit)          %u %u %u\n",
    sizeof (Unit), sizeof (Unit) / sizeof (int), UNITS (Unit, 1) )) ;
    DEBUG2 (("sizeof (Entry)         %u %u %u\n",
    sizeof (Entry), sizeof (Entry) / sizeof (int), UNITS (Entry, 1) )) ;
    DEBUG2 (("sizeof (Tuple)         %u %u %u\n",
    sizeof (Tuple), sizeof (Tuple) / sizeof (int), UNITS (Tuple, 1) )) ;
    DEBUG2 (("sizeof (Tuple *)       %u %u %u\n",
    sizeof (Tuple *), sizeof (Tuple *) / sizeof (int), UNITS (Tuple *, 1) )) ;
    DEBUG2 (("sizeof (Element)       %u %u %u\n",
    sizeof (Element), sizeof (Element) / sizeof (int), UNITS (Element, 1) )) ;
    DEBUG2 (("sizeof (Element *)     %u %u %u\n",
    sizeof (Element *), sizeof (Element *) / sizeof (int),
    UNITS (Element *, 1) )) ;
    DEBUG2 (("sizeof (WorkType)      %u %u %u\n",
    sizeof (WorkType), sizeof (WorkType) / sizeof (int),
    UNITS (WorkType, 1) )) ;
    DEBUG2 (("sizeof (NumericType)   %u %u %u\n",
    sizeof (NumericType), sizeof (NumericType) / sizeof (int),
    UNITS (NumericType, 1) )) ;
    DEBUG2 (("sizeof (SymbolicType)  %u %u %u\n",
    sizeof (SymbolicType), sizeof (SymbolicType) / sizeof (int),
    UNITS (SymbolicType, 1) )) ;

}


/* ========================================================================== */
/* === UMF_dump_rowmerge ==================================================== */
/* ========================================================================== */

GLOBAL void UMF_dump_rowmerge
(
    NumericType *Numeric,
    SymbolicType *Symbolic,
    WorkType *Work
)
{
    Int *Front_leftmostdesc, *Front_1strow, *Front_new1strow, row1, row2,
    fleftmost, nfr, n_row, *Row_degree, i, frontid, row ;

    nfr = Symbolic->nfr ;
    DEBUG3 (("\n================== Row merge sets: nfr "ID"\n", nfr)) ;
    Front_leftmostdesc = Symbolic->Front_leftmostdesc ;
    Front_1strow = Symbolic->Front_1strow ;
    Front_new1strow = Work->Front_new1strow ;
    n_row = Symbolic->n_row ;
    Row_degree = Numeric->Rperm ;
    frontid = Work->frontid ;

    for (i = frontid ; i <= nfr ; i++)
    {
    DEBUG3 (("----------------------\n")) ;
    if (i == nfr) DEBUG3 (("Dummy: ")) ;
    DEBUG3 (("Front "ID" 1strow "ID" new1strow "ID" leftmostdesc "ID,
        i, Front_1strow [i], Front_new1strow [i], Front_leftmostdesc [i])) ;
    DEBUG3 ((" parent "ID" pivcol "ID"\n", Symbolic->Front_parent [i],
        Symbolic->Front_npivcol [i])) ;

    if (i == nfr)
    {
        fleftmost = -1 ;
        row1 = Front_new1strow [i] ;
        row2 = n_row-1 ;
    }
    else
    {
        fleftmost = Front_leftmostdesc [i] ;
        row1 = Front_new1strow [fleftmost] ;
        row2 = Front_1strow [i+1] - 1 ;
    }
    DEBUG3 (("Leftmost: "ID"  Rows ["ID" to "ID"], search ["ID" to "ID"]\n",
        fleftmost, Front_1strow [i], row2, row1, row2)) ;

    for (row = row1 ; row <= row2 ; row++)
    {
        ASSERT (row >= 0 && row < n_row) ;
        DEBUG3 (("   Row "ID" live: %d\n", row, NON_PIVOTAL_ROW (row))) ;
    }
    }
}

/* ========================================================================== */
/* === UMF_dump_diagonal_map ================================================ */
/* ========================================================================== */

GLOBAL void UMF_dump_diagonal_map
(
    Int Diagonal_map [ ],
    Int Diagonal_imap [ ],
    Int n1,
    Int nn,
    Int nempty
)
{
    Int row, col ;
    if (Diagonal_map != (Int *) NULL)
    {
    DEBUG2 (("\nDump the Diagonal_map: n1 "ID" nn "ID" nempty "ID"\n",
        n1, nn, nempty)) ;
    for (col = n1 ; col < nn - nempty ; col++)
    {
        row = Diagonal_map [col] ;
        DEBUG2 (("     Diagonal_map [col = "ID"] gives "ID": ",
        col, row)) ;
        row = UNFLIP (row) ;
        DEBUG2 ((" row "ID"\n", row)) ;
        ASSERT (Diagonal_imap [row] == col) ;
    }
    }
}

#endif /* NDEBUG */