--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/liboctave/UMFPACK/UMFPACK/Source/umf_scale_column.c	Fri Feb 25 19:55:28 2005 +0000
@@ -0,0 +1,433 @@
+/* ========================================================================== */
+/* === UMF_scale_column ===================================================== */
+/* ========================================================================== */
+
+/* -------------------------------------------------------------------------- */
+/* 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                       */
+/* -------------------------------------------------------------------------- */
+
+/*
+    Scale the current pivot column, move the pivot row and column into place,
+    and log the permutation.
+*/
+
+#include "umf_internal.h"
+#include "umf_mem_free_tail_block.h"
+#include "umf_scale.h"
+
+/* ========================================================================== */
+/* === shift_pivot_row ====================================================== */
+/* ========================================================================== */
+
+/* Except for the BLAS, most of the time is typically spent in the following
+ * shift_pivot_row routine.  It copies the pivot row into the U block, and
+ * then fills in the whole in the C block by shifting the last row of C into
+ * the row vacated by the pivot row.
+ */
+
+PRIVATE void shift_pivot_row (Entry *Fd, Entry *Fs, Entry *Fe, Int len, Int d)
+{
+    Int j ;
+#pragma ivdep
+    for (j = 0 ; j < len ; j++)
+    {
+	Fd [j]   = Fs [j*d] ;
+	Fs [j*d] = Fe [j*d] ;
+    }
+}
+
+/* ========================================================================== */
+/* === UMF_scale_column ===================================================== */
+/* ========================================================================== */
+
+GLOBAL void UMF_scale_column
+(
+    NumericType *Numeric,
+    WorkType *Work
+)
+{
+    /* ---------------------------------------------------------------------- */
+    /* local variables */
+    /* ---------------------------------------------------------------------- */
+
+    Entry pivot_value ;
+    Entry *Fcol, *Flublock, *Flblock, *Fublock, *Fcblock ;
+    Int k, k1, fnr_curr, fnrows, fncols, *Frpos, *Fcpos, pivrow, pivcol,
+	*Frows, *Fcols, fnc_curr, fnpiv, *Row_tuples, nb,
+	*Col_tuples, *Rperm, *Cperm, fspos, col2, row2 ;
+#ifndef NDEBUG
+    Int *Col_degree, *Row_degree ;
+#endif
+
+    /* ---------------------------------------------------------------------- */
+    /* get parameters */
+    /* ---------------------------------------------------------------------- */
+
+    fnrows = Work->fnrows ;
+    fncols = Work->fncols ;
+    fnpiv = Work->fnpiv ;
+
+    /* ---------------------------------------------------------------------- */
+
+    Rperm = Numeric->Rperm ;
+    Cperm = Numeric->Cperm ;
+
+    /* ---------------------------------------------------------------------- */
+
+    Flublock = Work->Flublock ;
+    Flblock  = Work->Flblock ;
+    Fublock  = Work->Fublock ;
+    Fcblock  = Work->Fcblock ;
+
+    fnr_curr = Work->fnr_curr ;
+    fnc_curr = Work->fnc_curr ;
+    Frpos = Work->Frpos ;
+    Fcpos = Work->Fcpos ;
+    Frows = Work->Frows ;
+    Fcols = Work->Fcols ;
+    pivrow = Work->pivrow ;
+    pivcol = Work->pivcol ;
+
+    ASSERT (pivrow >= 0 && pivrow < Work->n_row) ;
+    ASSERT (pivcol >= 0 && pivcol < Work->n_col) ;
+
+#ifndef NDEBUG
+    Col_degree = Numeric->Cperm ;	/* for NON_PIVOTAL_COL macro */
+    Row_degree = Numeric->Rperm ;	/* for NON_PIVOTAL_ROW macro */
+#endif
+
+    Row_tuples = Numeric->Uip ;
+    Col_tuples = Numeric->Lip ;
+    nb = Work->nb ;
+
+#ifndef NDEBUG
+    ASSERT (fnrows == Work->fnrows_new + 1) ;
+    ASSERT (fncols == Work->fncols_new + 1) ;
+    DEBUG1 (("SCALE COL: fnrows "ID" fncols "ID"\n", fnrows, fncols)) ;
+    DEBUG2 (("\nFrontal matrix, including all space:\n"
+		"fnr_curr "ID" fnc_curr "ID" nb    "ID"\n"
+		"fnrows   "ID" fncols   "ID" fnpiv "ID"\n",
+		fnr_curr, fnc_curr, nb, fnrows, fncols, fnpiv)) ;
+    DEBUG2 (("\nJust the active part:\n")) ;
+    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, nb, fnpiv, fnpiv) ;
+#endif
+
+    /* ====================================================================== */
+    /* === Shift pivot row and column ======================================= */
+    /* ====================================================================== */
+
+    /* ---------------------------------------------------------------------- */
+    /* move pivot column into place */
+    /* ---------------------------------------------------------------------- */
+
+    /* Note that the pivot column is already in place.  Just shift the last
+     * column into the position vacated by the pivot column. */
+
+    fspos = Fcpos [pivcol] ;
+
+    /* one less column in the contribution block */
+    fncols = --(Work->fncols) ;
+
+    if (fspos != fncols * fnr_curr)
+    {
+
+	Int fs = fspos / fnr_curr ;
+
+	DEBUG6 (("Shift pivot column in front\n")) ;
+	DEBUG6 (("fspos: "ID" flpos: "ID"\n", fspos, fncols * fnr_curr)) ;
+
+	/* ------------------------------------------------------------------ */
+	/* move Fe => Fs */
+	/* ------------------------------------------------------------------ */
+
+	/* column of the contribution block: */
+	{
+	    /* Fs: current position of pivot column in contribution block */
+	    /* Fe: position of last column in contribution block */
+	    Int i ;
+	    Entry *Fs, *Fe ;
+	    Fs = Fcblock + fspos ;
+	    Fe = Fcblock + fncols * fnr_curr ;
+#pragma ivdep
+	    for (i = 0 ; i < fnrows ; i++)
+	    {
+		Fs [i] = Fe [i] ;
+	    }
+	}
+
+	/* column of the U2 block */
+	{
+	    /* Fs: current position of pivot column in U block */
+	    /* Fe: last column in U block */
+	    Int i ;
+	    Entry *Fs, *Fe ;
+	    Fs = Fublock + fs ;
+	    Fe = Fublock + fncols ;
+#pragma ivdep
+	    for (i = 0 ; i < fnpiv ; i++)
+	    {
+		Fs [i * fnc_curr] = Fe [i * fnc_curr] ;
+	    }
+	}
+
+	/* move column Fe to Fs in the Fcols pattern */
+	col2 = Fcols [fncols] ;
+	Fcols [fs] = col2 ;
+	Fcpos [col2] = fspos ;
+    }
+
+    /* pivot column is no longer in the frontal matrix */
+    Fcpos [pivcol] = EMPTY ;
+
+#ifndef NDEBUG
+    DEBUG2 (("\nFrontal matrix after col swap, including all space:\n"
+		"fnr_curr "ID" fnc_curr "ID" nb    "ID"\n"
+		"fnrows   "ID" fncols   "ID" fnpiv "ID"\n",
+		fnr_curr, fnc_curr, nb,
+		fnrows, fncols, fnpiv)) ;
+    DEBUG2 (("\nJust the active part:\n")) ;
+    DEBUG7 (("C  block: ")) ;
+    UMF_dump_dense (Fcblock,  fnr_curr, fnrows, fncols) ;
+    DEBUG7 (("L  block: ")) ;
+    UMF_dump_dense (Flblock,  fnr_curr, fnrows, fnpiv+1);
+    DEBUG7 (("U' block: ")) ;
+    UMF_dump_dense (Fublock,  fnc_curr, fncols, fnpiv) ;
+    DEBUG7 (("LU block: ")) ;
+    UMF_dump_dense (Flublock, nb, fnpiv, fnpiv+1) ;
+#endif
+
+    /* ---------------------------------------------------------------------- */
+    /* move pivot row into place */
+    /* ---------------------------------------------------------------------- */
+
+    fspos = Frpos [pivrow] ;
+
+    /* one less row in the contribution block */
+    fnrows = --(Work->fnrows) ;
+
+    DEBUG6 (("Swap/shift pivot row in front:\n")) ;
+    DEBUG6 (("fspos: "ID" flpos: "ID"\n", fspos, fnrows)) ;
+
+    if (fspos == fnrows)
+    {
+
+	/* ------------------------------------------------------------------ */
+	/* move Fs => Fd */
+	/* ------------------------------------------------------------------ */
+
+	DEBUG6 (("row case 1\n")) ;
+
+	/* row of the contribution block: */
+	{
+	    Int j ;
+	    Entry *Fd, *Fs ;
+	    Fd = Fublock + fnpiv * fnc_curr ;
+	    Fs = Fcblock + fspos ;
+#pragma ivdep
+	    for (j = 0 ; j < fncols ; j++)
+	    {
+		Fd [j] = Fs [j * fnr_curr] ;
+	    }
+	}
+
+	/* row of the L2 block: */
+	if (Work->pivrow_in_front)
+	{
+	    Int j ;
+	    Entry *Fd, *Fs ;
+	    Fd = Flublock + fnpiv ;
+	    Fs = Flblock  + fspos ;
+#pragma ivdep
+	    for (j = 0 ; j <= fnpiv ; j++)
+	    {
+		Fd [j * nb] = Fs [j * fnr_curr] ;
+	    }
+	}
+	else
+	{
+	    Int j ;
+	    Entry *Fd, *Fs ;
+	    Fd = Flublock + fnpiv ;
+	    Fs = Flblock  + fspos ;
+#pragma ivdep
+	    for (j = 0 ; j < fnpiv ; j++)
+	    {
+		ASSERT (IS_ZERO (Fs [j * fnr_curr])) ;
+		CLEAR (Fd [j * nb]) ;
+	    }
+	    Fd [fnpiv * nb] = Fs [fnpiv * fnr_curr] ;
+	}
+    }
+    else
+    {
+
+	/* ------------------------------------------------------------------ */
+	/* move Fs => Fd */
+	/* move Fe => Fs */
+	/* ------------------------------------------------------------------ */
+
+	DEBUG6 (("row case 2\n")) ;
+	/* this is the most common case, by far */
+
+	/* row of the contribution block: */
+	{
+	    /* Fd: destination of pivot row on U block */
+	    /* Fs: current position of pivot row in contribution block */
+	    /* Fe: position of last row in contribution block */
+	    Entry *Fd, *Fs, *Fe ;
+	    Fd = Fublock + fnpiv * fnc_curr ;
+	    Fs = Fcblock + fspos ;
+	    Fe = Fcblock + fnrows ;
+	    shift_pivot_row (Fd, Fs, Fe, fncols, fnr_curr) ;
+	}
+
+	/* row of the L2 block: */
+	if (Work->pivrow_in_front)
+	{
+	    /* Fd: destination of pivot row in LU block */
+	    /* Fs: current position of pivot row in L block */
+	    /* Fe: last row in L block */
+	    Int j ;
+	    Entry *Fd, *Fs, *Fe ;
+	    Fd = Flublock + fnpiv ;
+	    Fs = Flblock  + fspos ;
+	    Fe = Flblock  + fnrows ;
+#pragma ivdep
+	    for (j = 0 ; j <= fnpiv ; j++)
+	    {
+		Fd [j * nb]       = Fs [j * fnr_curr] ;
+		Fs [j * fnr_curr] = Fe [j * fnr_curr] ;
+	    }
+	}
+	else
+	{
+	    Int j ;
+	    Entry *Fd, *Fs, *Fe ;
+	    Fd = Flublock + fnpiv ;
+	    Fs = Flblock  + fspos ;
+	    Fe = Flblock  + fnrows ;
+#pragma ivdep
+	    for (j = 0 ; j < fnpiv ; j++)
+	    {
+		ASSERT (IS_ZERO (Fs [j * fnr_curr])) ;
+		CLEAR (Fd [j * nb]) ;
+		Fs [j * fnr_curr] = Fe [j * fnr_curr] ;
+	    }
+	    Fd [fnpiv * nb]       = Fs [fnpiv * fnr_curr] ;
+	    Fs [fnpiv * fnr_curr] = Fe [fnpiv * fnr_curr] ;
+	}
+
+	/* move row Fe to Fs in the Frows pattern */
+	row2 = Frows [fnrows] ;
+	Frows [fspos] = row2 ;
+	Frpos [row2] = fspos ;
+
+    }
+    /* pivot row is no longer in the frontal matrix */
+    Frpos [pivrow] = EMPTY ;
+
+#ifndef NDEBUG
+    DEBUG2 (("\nFrontal matrix after row swap, including all space:\n"
+		"fnr_curr "ID" fnc_curr "ID" nb    "ID"\n"
+		"fnrows   "ID" fncols   "ID" fnpiv "ID"\n",
+		Work->fnr_curr, Work->fnc_curr, Work->nb,
+		Work->fnrows, Work->fncols, Work->fnpiv)) ;
+    DEBUG2 (("\nJust the active part:\n")) ;
+    DEBUG7 (("C  block: ")) ;
+    UMF_dump_dense (Fcblock,  fnr_curr, fnrows, fncols) ;
+    DEBUG7 (("L  block: ")) ;
+    UMF_dump_dense (Flblock,  fnr_curr, fnrows, fnpiv+1);
+    DEBUG7 (("U' block: ")) ;
+    UMF_dump_dense (Fublock,  fnc_curr, fncols, fnpiv+1) ;
+    DEBUG7 (("LU block: ")) ;
+    UMF_dump_dense (Flublock, nb, fnpiv+1, fnpiv+1) ;
+#endif
+
+    /* ---------------------------------------------------------------------- */
+    /* Frpos [row] >= 0 for each row in pivot column pattern.   */
+    /* offset into pattern is given by:				*/
+    /* Frpos [row] == offset - 1				*/
+    /* Frpos [pivrow] is EMPTY */
+
+    /* Fcpos [col] >= 0 for each col in pivot row pattern.	*/
+    /* Fcpos [col] == (offset - 1) * fnr_curr			*/
+    /* Fcpos [pivcol] is EMPTY */
+
+    /* Fcols [0..fncols-1] is the pivot row pattern (excl pivot cols) */
+    /* Frows [0..fnrows-1] is the pivot col pattern (excl pivot rows) */
+
+    /* ====================================================================== */
+    /* === scale pivot column =============================================== */
+    /* ====================================================================== */
+
+    /* pivot column (except for pivot entry itself) */
+    Fcol = Flblock + fnpiv * fnr_curr ;
+    /* fnpiv-th pivot in frontal matrix located in Flublock (fnpiv, fnpiv) */
+    pivot_value = Flublock [fnpiv + fnpiv * nb] ;
+
+    /* this is the kth global pivot */
+    k = Work->npiv + fnpiv ;
+
+    DEBUG4 (("Pivot value: ")) ;
+    EDEBUG4 (pivot_value) ;
+    DEBUG4 (("\n")) ;
+
+    UMF_scale (fnrows, pivot_value, Fcol) ;
+
+    /* ---------------------------------------------------------------------- */
+    /* deallocate the pivot row and pivot column tuples */
+    /* ---------------------------------------------------------------------- */
+
+    UMF_mem_free_tail_block (Numeric, Row_tuples [pivrow]) ;
+    UMF_mem_free_tail_block (Numeric, Col_tuples [pivcol]) ;
+
+    Row_tuples [pivrow] = 0 ;
+    Col_tuples [pivcol] = 0 ;
+
+    DEBUG5 (("number of pivots prior to this one: "ID"\n", k)) ;
+    ASSERT (NON_PIVOTAL_ROW (pivrow)) ;
+    ASSERT (NON_PIVOTAL_COL (pivcol)) ;
+
+    /* save row and column inverse permutation */
+    k1 = ONES_COMPLEMENT (k) ;
+    Rperm [pivrow] = k1 ;			/* aliased with Row_degree */
+    Cperm [pivcol] = k1 ;			/* aliased with Col_degree */
+
+    ASSERT (!NON_PIVOTAL_ROW (pivrow)) ;
+    ASSERT (!NON_PIVOTAL_COL (pivcol)) ;
+
+    /* ---------------------------------------------------------------------- */
+    /* Keep track of the pivot order.  This is the kth pivot row and column. */
+    /* ---------------------------------------------------------------------- */
+
+    /* keep track of pivot rows and columns in the LU, L, and U blocks */
+    ASSERT (fnpiv < MAXNB) ;
+    Work->Pivrow [fnpiv] = pivrow ;
+    Work->Pivcol [fnpiv] = pivcol ;
+
+    /* ====================================================================== */
+    /* === one step in the factorization is done ============================ */
+    /* ====================================================================== */
+
+    /* One more step is done, except for pending updates to the U and C blocks
+     * of this frontal matrix.  Those are saved up, and applied by
+     * UMF_blas3_update when enough pivots have accumulated.   Also, the
+     * LU factors for these pending pivots have not yet been stored. */
+
+    Work->fnpiv++ ;
+
+#ifndef NDEBUG
+    DEBUG7 (("Current frontal matrix: (after pivcol scale)\n")) ;
+    UMF_dump_current_front (Numeric, Work, TRUE) ;
+#endif
+
+}