--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/liboctave/UMFPACK/UMFPACK/Source/umf_create_element.c	Fri Feb 25 19:55:28 2005 +0000
@@ -0,0 +1,604 @@
+/* ========================================================================== */
+/* === UMF_create_element =================================================== */
+/* ========================================================================== */
+
+/* -------------------------------------------------------------------------- */
+/* UMFPACK Version 4.4, 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                       */
+/* -------------------------------------------------------------------------- */
+
+/*
+    Factorization of a frontal matrix is complete.  Create a new element for
+    later assembly into a subsequent frontal matrix.  Returns TRUE if
+    successful, FALSE if out of memory.
+*/
+
+#include "umf_internal.h"
+#include "umf_mem_alloc_element.h"
+#include "umf_mem_alloc_tail_block.h"
+#include "umf_mem_free_tail_block.h"
+#include "umf_get_memory.h"
+
+/* ========================================================================== */
+/* === copy_column ========================================================== */
+/* ========================================================================== */
+
+PRIVATE void copy_column (Int len, Entry *X, Entry *Y)
+{
+    Int i ;
+#pragma ivdep
+    for (i = 0 ; i < len ; i++)
+    {
+	Y [i] = X [i] ;
+    }
+}
+
+/* ========================================================================== */
+/* === UMF_create_element =================================================== */
+/* ========================================================================== */
+
+GLOBAL Int UMF_create_element
+(
+    NumericType *Numeric,
+    WorkType *Work,
+    SymbolicType *Symbolic
+)
+{
+    /* ---------------------------------------------------------------------- */
+    /* local variables */
+    /* ---------------------------------------------------------------------- */
+
+    Int j, col, row, *Fcols, *Frows, fnrows, fncols, *Cols, len, needunits, t1,
+	t2, size, e, i, *E, *Fcpos, *Frpos, *Rows, eloc, fnr_curr, f,
+	got_memory, *Row_tuples, *Row_degree, *Row_tlen, *Col_tuples, max_mark,
+	*Col_degree, *Col_tlen, nn, n_row, n_col, r2, c2, do_Fcpos ;
+    Entry *C, *Fcol ;
+    Element *ep ;
+    Unit *p, *Memory ;
+    Tuple *tp, *tp1, *tp2, tuple, *tpend ;
+#ifndef NDEBUG
+    DEBUG2 (("FRONTAL WRAPUP\n")) ;
+    UMF_dump_current_front (Numeric, Work, TRUE) ;
+#endif
+
+    /* ---------------------------------------------------------------------- */
+    /* get parameters */
+    /* ---------------------------------------------------------------------- */
+
+    ASSERT (Work->fnpiv == 0) ;
+    ASSERT (Work->fnzeros == 0) ;
+    Row_degree = Numeric->Rperm ;
+    Row_tuples = Numeric->Uip ;
+    Row_tlen   = Numeric->Uilen ;
+    Col_degree = Numeric->Cperm ;
+    Col_tuples = Numeric->Lip ;
+    Col_tlen   = Numeric->Lilen ;
+    n_row = Work->n_row ;
+    n_col = Work->n_col ;
+    nn = MAX (n_row, n_col) ;
+    Fcols = Work->Fcols ;
+    Frows = Work->Frows ;
+    Fcpos = Work->Fcpos ;
+    Frpos = Work->Frpos ;
+    Memory = Numeric->Memory ;
+    fncols = Work->fncols ;
+    fnrows = Work->fnrows ;
+
+    tp = (Tuple *) NULL ;
+    tp1 = (Tuple *) NULL ;
+    tp2 = (Tuple *) NULL ;
+
+    /* ---------------------------------------------------------------------- */
+    /* add the current frontal matrix to the degrees of each column */
+    /* ---------------------------------------------------------------------- */
+
+    if (!Symbolic->fixQ)
+    {
+	/* but only if the column ordering is not fixed */
+#pragma ivdep
+	for (j = 0 ; j < fncols ; j++)
+	{
+	    /* add the current frontal matrix to the degree */
+	    ASSERT (Fcols [j] >= 0 && Fcols [j] < n_col) ;
+	    Col_degree [Fcols [j]] += fnrows ;
+	}
+    }
+
+    /* ---------------------------------------------------------------------- */
+    /* add the current frontal matrix to the degrees of each row */
+    /* ---------------------------------------------------------------------- */
+
+#pragma ivdep
+    for (i = 0 ; i < fnrows ; i++)
+    {
+	/* add the current frontal matrix to the degree */
+	ASSERT (Frows [i] >= 0 && Frows [i] < n_row) ;
+	Row_degree [Frows [i]] += fncols ;
+    }
+
+    /* ---------------------------------------------------------------------- */
+    /* Reset the external degree counters */
+    /* ---------------------------------------------------------------------- */
+
+    E = Work->E ;
+    max_mark = MAX_MARK (nn) ;
+
+    if (!Work->pivcol_in_front)
+    {
+	/* clear the external column degrees. no more Usons of current front */
+	Work->cdeg0 += (nn + 1) ;
+	if (Work->cdeg0 >= max_mark)
+	{
+	    /* guard against integer overflow.  This is very rare */
+	    DEBUG1 (("Integer overflow, cdeg\n")) ;
+	    Work->cdeg0 = 1 ;
+#pragma ivdep
+	    for (e = 1 ; e <= Work->nel ; e++)
+	    {
+		if (E [e])
+		{
+		    ep = (Element *) (Memory + E [e]) ;
+		    ep->cdeg = 0 ;
+		}
+	    }
+	}
+    }
+
+    if (!Work->pivrow_in_front)
+    {
+	/* clear the external row degrees.  no more Lsons of current front */
+	Work->rdeg0 += (nn + 1) ;
+	if (Work->rdeg0 >= max_mark)
+	{
+	    /* guard against integer overflow.  This is very rare */
+	    DEBUG1 (("Integer overflow, rdeg\n")) ;
+	    Work->rdeg0 = 1 ;
+#pragma ivdep
+	    for (e = 1 ; e <= Work->nel ; e++)
+	    {
+		if (E [e])
+		{
+		    ep = (Element *) (Memory + E [e]) ;
+		    ep->rdeg = 0 ;
+		}
+	    }
+	}
+    }
+
+    /* ---------------------------------------------------------------------- */
+    /* clear row/col offsets */
+    /* ---------------------------------------------------------------------- */
+
+    if (!Work->pivrow_in_front)
+    {
+#pragma ivdep
+	for (j = 0 ; j < fncols ; j++)
+	{
+	    Fcpos [Fcols [j]] = EMPTY ;
+	}
+    }
+
+    if (!Work->pivcol_in_front)
+    {
+#pragma ivdep
+	for (i = 0 ; i < fnrows ; i++)
+	{
+	    Frpos [Frows [i]] = EMPTY ;
+	}
+    }
+
+    if (fncols <= 0 || fnrows <= 0)
+    {
+	/* no element to create */
+	DEBUG2 (("Element evaporation\n")) ;
+	Work->prior_element = EMPTY ;
+	return (TRUE) ;
+    }
+
+    /* ---------------------------------------------------------------------- */
+    /* create element for later assembly */
+    /* ---------------------------------------------------------------------- */
+
+#ifndef NDEBUG
+    UMF_allocfail = FALSE ;
+    if (UMF_gprob > 0)
+    {
+	double rrr = ((double) (rand ( ))) / (((double) RAND_MAX) + 1) ;
+	DEBUG4 (("Check random %e %e\n", rrr, UMF_gprob)) ;
+	UMF_allocfail = rrr < UMF_gprob ;
+	if (UMF_allocfail) DEBUGm2 (("Random garbage collection (create)\n"));
+    }
+#endif
+
+    needunits = 0 ;
+    got_memory = FALSE ;
+    eloc = UMF_mem_alloc_element (Numeric, fnrows, fncols, &Rows, &Cols, &C,
+	&needunits, &ep) ;
+
+    /* if UMF_get_memory needs to be called */
+    if (Work->do_grow)
+    {
+	/* full compaction of current frontal matrix, since UMF_grow_front will
+	 * be called next anyway. */
+	r2 = fnrows ;
+	c2 = fncols ;
+	do_Fcpos = FALSE ;
+    }
+    else
+    {
+	/* partial compaction. */
+	r2 = MAX (fnrows, Work->fnrows_new + 1) ;
+	c2 = MAX (fncols, Work->fncols_new + 1) ;
+	/* recompute Fcpos if pivot row is in the front */
+	do_Fcpos = Work->pivrow_in_front ;
+    }
+
+    if (!eloc)
+    {
+	/* Do garbage collection, realloc, and try again. */
+	/* Compact the current front if it needs to grow anyway. */
+	/* Note that there are no pivot rows or columns in the current front */
+	DEBUGm3 (("get_memory from umf_create_element, 1\n")) ;
+	if (!UMF_get_memory (Numeric, Work, needunits, r2, c2, do_Fcpos))
+	{
+	    /* :: out of memory in umf_create_element (1) :: */
+	    DEBUGm4 (("out of memory: create element (1)\n")) ;
+	    return (FALSE) ;	/* out of memory */
+	}
+	got_memory = TRUE ;
+	Memory = Numeric->Memory ;
+	eloc = UMF_mem_alloc_element (Numeric, fnrows, fncols, &Rows, &Cols, &C,
+	    &needunits, &ep) ;
+	ASSERT (eloc >= 0) ;
+	if (!eloc)
+	{
+	    /* :: out of memory in umf_create_element (2) :: */
+	    DEBUGm4 (("out of memory: create element (2)\n")) ;
+	    return (FALSE) ;	/* out of memory */
+	}
+    }
+
+    e = ++(Work->nel) ;	/* get the name of this new frontal matrix */
+    Work->prior_element = e ;
+    DEBUG8 (("wrapup e "ID" nel "ID"\n", e, Work->nel)) ;
+
+    ASSERT (e > 0 && e < Work->elen) ;
+    ASSERT (E [e] == 0) ;
+    E [e] = eloc ;
+
+    if (Work->pivcol_in_front)
+    {
+	/* the new element is a Uson of the next frontal matrix */
+	ep->cdeg = Work->cdeg0 ;
+    }
+
+    if (Work->pivrow_in_front)
+    {
+	/* the new element is an Lson of the next frontal matrix */
+	ep->rdeg = Work->rdeg0 ;
+    }
+
+    /* ---------------------------------------------------------------------- */
+    /* copy frontal matrix into the new element */
+    /* ---------------------------------------------------------------------- */
+
+#pragma ivdep
+    for (i = 0 ; i < fnrows ; i++)
+    {
+	Rows [i] = Frows [i] ;
+    }
+#pragma ivdep
+    for (i = 0 ; i < fncols ; i++)
+    {
+	Cols [i] = Fcols [i] ;
+    }
+    Fcol = Work->Fcblock ;
+    DEBUG0 (("copy front "ID" by "ID"\n", fnrows, fncols)) ;
+    fnr_curr = Work->fnr_curr ;
+    ASSERT (fnr_curr >= 0 && fnr_curr % 2 == 1) ;
+    for (j = 0 ; j < fncols ; j++)
+    {
+	copy_column (fnrows, Fcol, C) ;
+#if 0
+#ifdef USE_NO_BLAS
+	copy_column (fnrows, Fcol, C) ;
+#else
+	could also use BLAS-COPY (fnrows, Fcol, C) here, but it is typically
+	not as fast as the inlined copy_column subroutine, above.
+#endif
+	for (i = 0 ; i < fnrows ; i++)
+	{
+	    C [i] = Fcol [i] ;
+	}
+#endif
+	Fcol += fnr_curr ;
+	C += fnrows ;
+    }
+
+    DEBUG8 (("element copied\n")) ;
+
+    /* ---------------------------------------------------------------------- */
+    /* add tuples for the new element */
+    /* ---------------------------------------------------------------------- */
+
+    tuple.e = e ;
+
+    if (got_memory)
+    {
+
+	/* ------------------------------------------------------------------ */
+	/* UMF_get_memory ensures enough space exists for each new tuple */
+	/* ------------------------------------------------------------------ */
+
+	/* place (e,f) in the element list of each column */
+	for (tuple.f = 0 ; tuple.f < fncols ; tuple.f++)
+	{
+	    col = Fcols [tuple.f] ;
+	    ASSERT (col >= 0 && col < n_col) ;
+	    ASSERT (NON_PIVOTAL_COL (col)) ;
+	    ASSERT (Col_tuples [col]) ;
+	    tp = ((Tuple *) (Memory + Col_tuples [col])) + Col_tlen [col]++ ;
+	    *tp = tuple ;
+	}
+
+	/* ------------------------------------------------------------------ */
+
+	/* place (e,f) in the element list of each row */
+	for (tuple.f = 0 ; tuple.f < fnrows ; tuple.f++)
+	{
+	    row = Frows [tuple.f] ;
+	    ASSERT (row >= 0 && row < n_row) ;
+	    ASSERT (NON_PIVOTAL_ROW (row)) ;
+	    ASSERT (Row_tuples [row]) ;
+	    tp = ((Tuple *) (Memory + Row_tuples [row])) + Row_tlen [row]++ ;
+	    *tp = tuple ;
+	}
+
+    }
+    else
+    {
+
+	/* ------------------------------------------------------------------ */
+	/* place (e,f) in the element list of each column */
+	/* ------------------------------------------------------------------ */
+
+	/* might not have enough space for each tuple */
+
+	for (tuple.f = 0 ; tuple.f < fncols ; tuple.f++)
+	{
+	    col = Fcols [tuple.f] ;
+	    ASSERT (col >= 0 && col < n_col) ;
+	    ASSERT (NON_PIVOTAL_COL (col)) ;
+	    t1 = Col_tuples [col] ;
+	    DEBUG1 (("Placing on col:"ID" , tuples at "ID"\n",
+		col, Col_tuples [col])) ;
+
+	    size = 0 ;
+	    len = 0 ;
+
+	    if (t1)
+	    {
+		p = Memory + t1 ;
+		tp = (Tuple *) p ;
+		size = GET_BLOCK_SIZE (p) ;
+		len = Col_tlen [col] ;
+		tp2 = tp + len ;
+	    }
+
+	    needunits = UNITS (Tuple, len + 1) ;
+	    DEBUG1 (("len: "ID" size: "ID" needunits: "ID"\n",
+		len, size, needunits));
+
+	    if (needunits > size && t1)
+	    {
+		/* prune the tuples */
+		tp1 = tp ;
+		tp2 = tp ;
+		tpend = tp + len ;
+		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 ;	/* already assembled */
+		    ASSERT (col == Cols [f]) ;
+		    *tp2++ = *tp ;	/* leave the tuple in the list */
+		}
+		len = tp2 - tp1 ;
+		Col_tlen [col] = len ;
+		needunits = UNITS (Tuple, len + 1) ;
+	    }
+
+	    if (needunits > size)
+	    {
+		/* no room exists - reallocate elsewhere */
+		DEBUG1 (("REALLOCATE Col: "ID", size "ID" to "ID"\n",
+		    col, size, 2*needunits)) ;
+
+#ifndef NDEBUG
+		UMF_allocfail = FALSE ;
+		if (UMF_gprob > 0)  /* a double relop, but ignore NaN case */
+		{
+		    double rrr = ((double) (rand ( ))) /
+			(((double) RAND_MAX) + 1) ;
+		    DEBUG1 (("Check random %e %e\n", rrr, UMF_gprob)) ;
+		    UMF_allocfail = rrr < UMF_gprob ;
+		    if (UMF_allocfail) DEBUGm2 (("Random gar. (col tuple)\n")) ;
+		}
+#endif
+
+		needunits = MIN (2*needunits, (Int) UNITS (Tuple, nn)) ;
+		t2 = UMF_mem_alloc_tail_block (Numeric, needunits) ;
+		if (!t2)
+		{
+		    /* :: get memory in umf_create_element (1) :: */
+		    /* get memory, reconstruct all tuple lists, and return */
+		    /* Compact the current front if it needs to grow anyway. */
+		    /* Note: no pivot rows or columns in the current front */
+		    DEBUGm4 (("get_memory from umf_create_element, 1\n")) ;
+		    return (UMF_get_memory (Numeric, Work, 0, r2, c2,do_Fcpos));
+		}
+		Col_tuples [col] = t2 ;
+		tp2 = (Tuple *) (Memory + t2) ;
+		if (t1)
+		{
+		    for (i = 0 ; i < len ; i++)
+		    {
+			*tp2++ = *tp1++ ;
+		    }
+		    UMF_mem_free_tail_block (Numeric, t1) ;
+		}
+	    }
+
+	    /* place the new (e,f) tuple in the element list of the column */
+	    Col_tlen [col]++ ;
+	    *tp2 = tuple ;
+	}
+
+	/* ------------------------------------------------------------------ */
+	/* place (e,f) in the element list of each row */
+	/* ------------------------------------------------------------------ */
+
+	for (tuple.f = 0 ; tuple.f < fnrows ; tuple.f++)
+	{
+	    row = Frows [tuple.f] ;
+	    ASSERT (row >= 0 && row < n_row) ;
+	    ASSERT (NON_PIVOTAL_ROW (row)) ;
+	    t1 = Row_tuples [row] ;
+	    DEBUG1 (("Placing on row:"ID" , tuples at "ID"\n",
+		row, Row_tuples [row])) ;
+
+	    size = 0 ;
+	    len = 0 ;
+	    if (t1)
+	    {
+		p = Memory + t1 ;
+		tp = (Tuple *) p ;
+		size = GET_BLOCK_SIZE (p) ;
+		len = Row_tlen [row] ;
+		tp2 = tp + len ;
+	    }
+
+	    needunits = UNITS (Tuple, len + 1) ;
+	    DEBUG1 (("len: "ID" size: "ID" needunits: "ID"\n",
+		len, size, needunits)) ;
+
+	    if (needunits > size && t1)
+	    {
+		/* prune the tuples */
+		tp1 = tp ;
+		tp2 = tp ;
+		tpend = tp + len ;
+		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 ;	/* already assembled */
+		    ASSERT (row == Rows [f]) ;
+		    *tp2++ = *tp ;	/* leave the tuple in the list */
+		}
+		len = tp2 - tp1 ;
+		Row_tlen [row] = len ;
+		needunits = UNITS (Tuple, len + 1) ;
+	    }
+
+	    if (needunits > size)
+	    {
+		/* no room exists - reallocate elsewhere */
+		DEBUG1 (("REALLOCATE Row: "ID", size "ID" to "ID"\n",
+		    row, size, 2*needunits)) ;
+
+#ifndef NDEBUG
+		UMF_allocfail = FALSE ;
+		if (UMF_gprob > 0)  /* a double relop, but ignore NaN case */
+		{
+		    double rrr = ((double) (rand ( ))) /
+			(((double) RAND_MAX) + 1) ;
+		    DEBUG1 (("Check random %e %e\n", rrr, UMF_gprob)) ;
+		    UMF_allocfail = rrr < UMF_gprob ;
+		    if (UMF_allocfail) DEBUGm2 (("Random gar. (row tuple)\n")) ;
+		}
+#endif
+
+		needunits = MIN (2*needunits, (Int) UNITS (Tuple, nn)) ;
+		t2 = UMF_mem_alloc_tail_block (Numeric, needunits) ;
+		if (!t2)
+		{
+		    /* :: get memory in umf_create_element (2) :: */
+		    /* get memory, reconstruct all tuple lists, and return */
+		    /* Compact the current front if it needs to grow anyway. */
+		    /* Note: no pivot rows or columns in the current front */
+		    DEBUGm4 (("get_memory from umf_create_element, 2\n")) ;
+		    return (UMF_get_memory (Numeric, Work, 0, r2, c2,do_Fcpos));
+		}
+		Row_tuples [row] = t2 ;
+		tp2 = (Tuple *) (Memory + t2) ;
+		if (t1)
+		{
+		    for (i = 0 ; i < len ; i++)
+		    {
+			*tp2++ = *tp1++ ;
+		    }
+		    UMF_mem_free_tail_block (Numeric, t1) ;
+		}
+	    }
+
+	    /* place the new (e,f) tuple in the element list of the row */
+	    Row_tlen [row]++ ;
+	    *tp2 = tuple ;
+	}
+
+    }
+
+    /* ---------------------------------------------------------------------- */
+
+#ifndef NDEBUG
+    DEBUG1 (("Done extending\nFINAL: element row pattern: len="ID"\n", fncols));
+    for (j = 0 ; j < fncols ; j++) DEBUG1 ((""ID"\n", Fcols [j])) ;
+    DEBUG1 (("FINAL: element col pattern:  len="ID"\n", fnrows)) ;
+    for (j = 0 ; j < fnrows ; j++) DEBUG1 ((""ID"\n", Frows [j])) ;
+    for (j = 0 ; j < fncols ; j++)
+    {
+	col = Fcols [j] ;
+	ASSERT (col >= 0 && col < n_col) ;
+	UMF_dump_rowcol (1, Numeric, Work, col, !Symbolic->fixQ) ;
+    }
+    for (j = 0 ; j < fnrows ; j++)
+    {
+	row = Frows [j] ;
+	ASSERT (row >= 0 && row < n_row) ;
+	UMF_dump_rowcol (0, Numeric, Work, row, TRUE) ;
+    }
+    if (n_row < 1000 && n_col < 1000)
+    {
+	UMF_dump_memory (Numeric) ;
+    }
+    DEBUG1 (("New element, after filling with stuff: "ID"\n", e)) ;
+    UMF_dump_element (Numeric, Work, e, TRUE) ;
+    if (nn < 1000)
+    {
+	DEBUG4 (("Matrix dump, after New element: "ID"\n", e)) ;
+	UMF_dump_matrix (Numeric, Work, TRUE) ;
+    }
+    DEBUG3 (("FRONTAL WRAPUP DONE\n")) ;
+#endif
+
+    return (TRUE) ;
+}