--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/liboctave/UMFPACK/UMFPACK/Source/umf_grow_front.c	Fri Feb 25 19:55:28 2005 +0000
@@ -0,0 +1,293 @@
+/* ========================================================================== */
+/* === UMF_grow_front ======================================================= */
+/* ========================================================================== */
+
+/* -------------------------------------------------------------------------- */
+/* 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                       */
+/* -------------------------------------------------------------------------- */
+
+/* Current frontal matrix is too small.  Make it bigger. */
+
+#include "umf_internal.h"
+#include "umf_mem_free_tail_block.h"
+#include "umf_mem_alloc_tail_block.h"
+#include "umf_get_memory.h"
+
+GLOBAL Int UMF_grow_front
+(
+    NumericType *Numeric,
+    Int fnr2,		/* desired size is fnr2-by-fnc2 */
+    Int fnc2,
+    WorkType *Work,
+    Int do_what		/* -1: UMF_start_front
+			 * 0:  UMF_init_front, do not recompute Fcpos
+			 * 1:  UMF_extend_front
+			 * 2:  UMF_init_front, recompute Fcpos */
+)
+{
+    /* ---------------------------------------------------------------------- */
+    /* local variables */
+    /* ---------------------------------------------------------------------- */
+
+    double s ;
+    Entry *Fcold, *Fcnew ;
+    Int j, i, col, *Fcpos, *Fcols, fnrows_max, fncols_max, fnr_curr, nb,
+	fnrows_new, fncols_new, fnr_min, fnc_min, minsize,
+	newsize, fnrows, fncols, *E, eloc ;
+
+    /* ---------------------------------------------------------------------- */
+    /* get parameters */
+    /* ---------------------------------------------------------------------- */
+
+#ifndef NDEBUG
+    if (do_what != -1) UMF_debug++ ;
+    DEBUG0 (("\n\n====================GROW FRONT: do_what: "ID"\n", do_what)) ;
+    if (do_what != -1) UMF_debug-- ;
+    ASSERT (Work->do_grow) ;
+    ASSERT (Work->fnpiv == 0) ;
+#endif
+
+    Fcols = Work->Fcols ;
+    Fcpos = Work->Fcpos ;
+    E = Work->E ;
+
+    /* ---------------------------------------------------------------------- */
+    /* The current front is too small, find the new size */
+    /* ---------------------------------------------------------------------- */
+
+    /* maximum size of frontal matrix for this chain */
+    nb = Work->nb ;
+    fnrows_max = Work->fnrows_max + nb ;
+    fncols_max = Work->fncols_max + nb ;
+    ASSERT (fnrows_max >= 0 && (fnrows_max % 2) == 1) ;
+    DEBUG0 (("Max     size: "ID"-by-"ID" (incl. "ID" pivot block\n",
+	fnrows_max, fncols_max, nb)) ;
+
+    /* current dimensions of frontal matrix: fnr-by-fnc */
+    DEBUG0 (("Current : "ID"-by-"ID" (excl "ID" pivot blocks)\n",
+		Work->fnr_curr, Work->fnc_curr, nb)) ;
+    ASSERT (Work->fnr_curr >= 0) ;
+    ASSERT ((Work->fnr_curr % 2 == 1) || Work->fnr_curr == 0) ;
+
+    /* required dimensions of frontal matrix: fnr_min-by-fnc_min */
+    fnrows_new = Work->fnrows_new + 1 ;
+    fncols_new = Work->fncols_new + 1 ;
+    ASSERT (fnrows_new >= 0) ;
+    if (fnrows_new % 2 == 0) fnrows_new++ ;
+    fnrows_new += nb ;
+    fncols_new += nb ;
+    fnr_min = MIN (fnrows_new, fnrows_max) ;
+    fnc_min = MIN (fncols_new, fncols_max) ;
+    minsize = fnr_min * fnc_min ;
+    if (INT_OVERFLOW ((double) fnr_min * (double) fnc_min * sizeof (Entry)))
+    {
+	/* :: the minimum front size is bigger than the integer maximum :: */
+	return (FALSE) ;
+    }
+    ASSERT (fnr_min >= 0) ;
+    ASSERT (fnr_min % 2 == 1) ;
+
+    DEBUG0 (("Min     : "ID"-by-"ID"\n", fnr_min, fnc_min)) ;
+
+    /* grow the front to fnr2-by-fnc2, but no bigger than the maximum,
+     * and no smaller than the minumum. */
+    DEBUG0 (("Desired : ("ID"+"ID")-by-("ID"+"ID")\n", fnr2, nb, fnc2, nb)) ;
+    fnr2 += nb ;
+    fnc2 += nb ;
+    ASSERT (fnr2 >= 0) ;
+    if (fnr2 % 2 == 0) fnr2++ ;
+    fnr2 = MAX (fnr2, fnr_min) ;
+    fnc2 = MAX (fnc2, fnc_min) ;
+    fnr2 = MIN (fnr2, fnrows_max) ;
+    fnc2 = MIN (fnc2, fncols_max) ;
+    DEBUG0 (("Try     : "ID"-by-"ID"\n", fnr2, fnc2)) ;
+    ASSERT (fnr2 >= 0) ;
+    ASSERT (fnr2 % 2 == 1) ;
+
+    s = ((double) fnr2) * ((double) fnc2) ;
+    if (INT_OVERFLOW (s * sizeof (Entry)))
+    {
+	/* :: frontal matrix size int overflow :: */
+	/* the desired front size is bigger than the integer maximum */
+	/* compute a such that a*a*s < Int_MAX / sizeof (Entry) */
+	double a = 0.9 * sqrt ((Int_MAX / sizeof (Entry)) / s) ;
+	fnr2 = MAX (fnr_min, a * fnr2) ;
+	fnc2 = MAX (fnc_min, a * fnc2) ;
+	/* the new frontal size is a*r*a*c = a*a*s */
+	newsize = fnr2 * fnc2 ;
+	ASSERT (fnr2 >= 0) ;
+	if (fnr2 % 2 == 0) fnr2++ ;
+	fnc2 = newsize / fnr2 ;
+    }
+
+    fnr2 = MAX (fnr2, fnr_min) ;
+    fnc2 = MAX (fnc2, fnc_min) ;
+    newsize = fnr2 * fnc2 ;
+
+    ASSERT (fnr2 >= 0) ;
+    ASSERT (fnr2 % 2 == 1) ;
+    ASSERT (fnr2 >= fnr_min) ;
+    ASSERT (fnc2 >= fnc_min) ;
+    ASSERT (newsize >= minsize) ;
+
+    /* ---------------------------------------------------------------------- */
+    /* free the current front if it is empty of any numerical values */
+    /* ---------------------------------------------------------------------- */
+
+    if (E [0] && do_what != 1)
+    {
+	/* free the current front, if it exists and has nothing in it */
+	DEBUG0 (("Freeing empty front\n")) ;
+	UMF_mem_free_tail_block (Numeric, E [0]) ;
+	E [0] = 0 ;
+	Work->Flublock = (Entry *) NULL ;
+	Work->Flblock  = (Entry *) NULL ;
+	Work->Fublock  = (Entry *) NULL ;
+	Work->Fcblock  = (Entry *) NULL ;
+    }
+
+    /* ---------------------------------------------------------------------- */
+    /* allocate the new front, doing garbage collection if necessary */
+    /* ---------------------------------------------------------------------- */
+
+#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 garbage collection (grow)\n")) ;
+    }
+#endif
+
+    DEBUG0 (("Attempt size: "ID"-by-"ID"\n", fnr2, fnc2)) ;
+    eloc = UMF_mem_alloc_tail_block (Numeric, UNITS (Entry, newsize)) ;
+
+    if (!eloc)
+    {
+	/* Do garbage collection, realloc, and try again. Compact the current
+	 * contribution block in the front to fnrows-by-fncols.  Note that
+	 * there are no pivot rows/columns in current front.  Do not recompute
+	 * Fcpos in UMF_garbage_collection. */
+	DEBUGm3 (("get_memory from umf_grow_front\n")) ;
+	if (!UMF_get_memory (Numeric, Work, 1 + UNITS (Entry, newsize),
+	    Work->fnrows, Work->fncols, FALSE))
+	{
+	    /* :: out of memory in umf_grow_front :: */
+	    return (FALSE) ;	/* out of memory */
+	}
+	DEBUG0 (("Attempt size: "ID"-by-"ID" again\n", fnr2, fnc2)) ;
+	eloc = UMF_mem_alloc_tail_block (Numeric, UNITS (Entry, newsize)) ;
+    }
+
+    /* try again with something smaller */
+    while ((fnr2 != fnr_min || fnc2 != fnc_min) && !eloc)
+    {
+	fnr2 = MIN (fnr2 - 2, fnr2 * UMF_REALLOC_REDUCTION) ;
+	fnc2 = MIN (fnc2 - 2, fnc2 * UMF_REALLOC_REDUCTION) ;
+	ASSERT (fnr_min >= 0) ;
+	ASSERT (fnr_min % 2 == 1) ;
+	fnr2 = MAX (fnr_min, fnr2) ;
+	fnc2 = MAX (fnc_min, fnc2) ;
+	ASSERT (fnr2 >= 0) ;
+	if (fnr2 % 2 == 0) fnr2++ ;
+	newsize = fnr2 * fnc2 ;
+	DEBUGm3 (("Attempt smaller size: "ID"-by-"ID" minsize "ID"-by-"ID"\n",
+	    fnr2, fnc2, fnr_min, fnc_min)) ;
+	eloc = UMF_mem_alloc_tail_block (Numeric, UNITS (Entry, newsize)) ;
+    }
+
+    /* try again with the smallest possible size */
+    if (!eloc)
+    {
+	fnr2 = fnr_min ;
+	fnc2 = fnc_min ;
+	newsize = minsize ;
+	DEBUG0 (("Attempt minsize: "ID"-by-"ID"\n", fnr2, fnc2)) ;
+	eloc = UMF_mem_alloc_tail_block (Numeric, UNITS (Entry, newsize)) ;
+    }
+
+    if (!eloc)
+    {
+	/* out of memory */
+	return (FALSE) ;
+    }
+
+    ASSERT (fnr2 >= 0) ;
+    ASSERT (fnr2 % 2 == 1) ;
+    ASSERT (fnr2 >= fnr_min && fnc2 >= fnc_min) ;
+
+    /* ---------------------------------------------------------------------- */
+    /* copy the old frontal matrix into the new one */
+    /* ---------------------------------------------------------------------- */
+
+    /* old contribution block (if any) */
+    fnr_curr = Work->fnr_curr ;	    /* garbage collection can change fn*_curr */
+    ASSERT (fnr_curr >= 0) ;
+    ASSERT ((fnr_curr % 2 == 1) || fnr_curr == 0) ;
+    fnrows = Work->fnrows ;
+    fncols = Work->fncols ;
+    Fcold = Work->Fcblock ;
+
+    /* remove nb from the sizes */
+    fnr2 -= nb ;
+    fnc2 -= nb ;
+
+    /* new frontal matrix */
+    Work->Flublock = (Entry *) (Numeric->Memory + eloc) ;
+    Work->Flblock  = Work->Flublock + nb * nb ;
+    Work->Fublock  = Work->Flblock  + nb * fnr2 ;
+    Work->Fcblock  = Work->Fublock  + nb * fnc2 ;
+    Fcnew = Work->Fcblock ;
+
+    if (E [0])
+    {
+	/* copy the old contribution block into the new one */
+	for (j = 0 ; j < fncols ; j++)
+	{
+	    col = Fcols [j] ;
+	    DEBUG1 (("copy col "ID" \n",col)) ;
+	    ASSERT (col >= 0 && col < Work->n_col) ;
+	    for (i = 0 ; i < fnrows ; i++)
+	    {
+		Fcnew [i] = Fcold [i] ;
+	    }
+	    Fcnew += fnr2 ;
+	    Fcold += fnr_curr ;
+	    DEBUG1 (("new offset col "ID" "ID"\n",col, j * fnr2)) ;
+	    Fcpos [col] = j * fnr2 ;
+	}
+    }
+    else if (do_what == 2)
+    {
+	/* just find the new column offsets */
+	for (j = 0 ; j < fncols ; j++)
+	{
+	    col = Fcols [j] ;
+	    DEBUG1 (("new offset col "ID" "ID"\n",col, j * fnr2)) ;
+	    Fcpos [col] = j * fnr2 ;
+	}
+    }
+
+    /* free the old frontal matrix */
+    UMF_mem_free_tail_block (Numeric, E [0]) ;
+
+    /* ---------------------------------------------------------------------- */
+    /* new frontal matrix size */
+    /* ---------------------------------------------------------------------- */
+
+    E [0] = eloc ;
+    Work->fnr_curr = fnr2 ;	    /* C block is fnr2-by-fnc2 */
+    Work->fnc_curr = fnc2 ;
+    Work->fcurr_size = newsize ;    /* including LU, L, U, and C blocks */
+    Work->do_grow = FALSE ;	    /* the front has just been grown */
+
+    ASSERT (Work->fnr_curr >= 0) ;
+    ASSERT (Work->fnr_curr % 2 == 1) ;
+    DEBUG0 (("Newly grown front: "ID"+"ID" by "ID"+"ID"\n", Work->fnr_curr,
+	nb, Work->fnc_curr, nb)) ;
+    return (TRUE) ;
+}