--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/liboctave/UMFPACK/UMFPACK/Source/umf_kernel_wrapup.c	Fri Feb 25 19:55:28 2005 +0000
@@ -0,0 +1,466 @@
+/* ========================================================================== */
+/* === UMF_kernel_wrapup ==================================================== */
+/* ========================================================================== */
+
+/* -------------------------------------------------------------------------- */
+/* 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                       */
+/* -------------------------------------------------------------------------- */
+
+/* The matrix is factorized.  Finish the LU data structure. */
+
+#include "umf_internal.h"
+
+GLOBAL void UMF_kernel_wrapup
+(
+    NumericType *Numeric,
+    SymbolicType *Symbolic,
+    WorkType *Work
+)
+{
+
+    /* ---------------------------------------------------------------------- */
+    /* local variables */
+    /* ---------------------------------------------------------------------- */
+
+    Entry pivot_value ;
+    double d ;
+    Entry *D ;
+    Int i, k, col, row, llen, ulen, *ip, *Rperm, *Cperm, *Lilen, npiv, lp,
+	*Uilen, *Lip, *Uip, *Cperm_init, up, pivrow, pivcol, *Lpos, *Upos, *Wr,
+	*Wc, *Wp, *Frpos, *Fcpos, *Row_degree, *Col_degree, *Rperm_init,
+	n_row, n_col, n_inner, zero_pivot, nan_pivot, n1 ;
+
+#ifndef NDEBUG
+    UMF_dump_matrix (Numeric, Work, FALSE) ;
+#endif
+
+    DEBUG0 (("Kernel complete, Starting Kernel wrapup\n")) ;
+    n_row = Symbolic->n_row ;
+    n_col = Symbolic->n_col ;
+    n_inner = MIN (n_row, n_col) ;
+    Rperm = Numeric->Rperm ;
+    Cperm = Numeric->Cperm ;
+    Lilen = Numeric->Lilen ;
+    Uilen = Numeric->Uilen ;
+    Upos = Numeric->Upos ;
+    Lpos = Numeric->Lpos ;
+    Lip = Numeric->Lip ;
+    Uip = Numeric->Uip ;
+    D = Numeric->D ;
+
+    npiv = Work->npiv ;
+    Numeric->npiv = npiv ;
+    Numeric->ulen = Work->ulen ;
+
+    ASSERT (n_row == Numeric->n_row) ;
+    ASSERT (n_col == Symbolic->n_col) ;
+    DEBUG0 (("Wrap-up: npiv "ID" ulen "ID"\n", npiv, Numeric->ulen)) ;
+    ASSERT (npiv <= n_inner) ;
+
+    /* this will be nonzero only if matrix is singular or rectangular */
+    ASSERT (IMPLIES (npiv == n_col, Work->ulen == 0)) ;
+
+    /* ---------------------------------------------------------------------- */
+    /* find the smallest and largest entries in D */
+    /* ---------------------------------------------------------------------- */
+
+    for (k = 0 ; k < npiv ; k++)
+    {
+	pivot_value = D [k] ;
+	ABS (d, pivot_value) ;
+	zero_pivot = SCALAR_IS_ZERO (d) ;
+	nan_pivot = SCALAR_IS_NAN (d) ;
+
+	if (!zero_pivot)
+	{
+	    /* the pivot is nonzero, but might be Inf or NaN */
+	    Numeric->nnzpiv++ ;
+	}
+
+	if (k == 0)
+	{
+	    Numeric->min_udiag = d ;
+	    Numeric->max_udiag = d ;
+	}
+	else
+	{
+	    /* min (abs (diag (U))) behaves as follows:  If any entry is zero,
+	       then the result is zero (regardless of the presence of NaN's).
+	       Otherwise, if any entry is NaN, then the result is NaN.
+	       Otherwise, the result is the smallest absolute value on the
+	       diagonal of U.
+	    */
+
+	    if (SCALAR_IS_NONZERO (Numeric->min_udiag))
+	    {
+		if (zero_pivot || nan_pivot)
+		{
+		    Numeric->min_udiag = d ;
+		}
+		else if (!SCALAR_IS_NAN (Numeric->min_udiag))
+		{
+		    /* d and min_udiag are both non-NaN */
+		    Numeric->min_udiag = MIN (Numeric->min_udiag, d) ;
+		}
+	    }
+
+	    /*
+	       max (abs (diag (U))) behaves as follows:  If any entry is NaN
+	       then the result is NaN.  Otherise, the result is the largest
+	       absolute value on the diagonal of U.
+	    */
+
+	    if (nan_pivot)
+	    {
+		Numeric->max_udiag = d ;
+	    }
+	    else if (!SCALAR_IS_NAN (Numeric->max_udiag))
+	    {
+		/* d and max_udiag are both non-NaN */
+		Numeric->max_udiag = MAX (Numeric->max_udiag, d) ;
+	    }
+	}
+    }
+
+    /* ---------------------------------------------------------------------- */
+    /* check if matrix is singular or rectangular */
+    /* ---------------------------------------------------------------------- */
+
+    Col_degree = Cperm ;	/* for NON_PIVOTAL_COL macro */
+    Row_degree = Rperm ;	/* for NON_PIVOTAL_ROW macro */
+
+    if (npiv < n_row)
+    {
+	/* finalize the row permutation */
+	k = npiv ;
+	DEBUGm3 (("Singular pivot rows "ID" to "ID"\n", k, n_row-1)) ;
+	for (row = 0 ; row < n_row ; row++)
+	{
+	    if (NON_PIVOTAL_ROW (row))
+	    {
+		Rperm [row] = ONES_COMPLEMENT (k) ;
+		DEBUGm3 (("Singular row "ID" is k: "ID" pivot row\n", row, k)) ;
+		ASSERT (!NON_PIVOTAL_ROW (row)) ;
+		Lpos [row] = EMPTY ;
+		Uip [row] = EMPTY ;
+		Uilen [row] = 0 ;
+		k++ ;
+	    }
+	}
+	ASSERT (k == n_row) ;
+    }
+
+    if (npiv < n_col)
+    {
+	/* finalize the col permutation */
+	k = npiv ;
+	DEBUGm3 (("Singular pivot cols "ID" to "ID"\n", k, n_col-1)) ;
+	for (col = 0 ; col < n_col ; col++)
+	{
+	    if (NON_PIVOTAL_COL (col))
+	    {
+		Cperm [col] = ONES_COMPLEMENT (k) ;
+		DEBUGm3 (("Singular col "ID" is k: "ID" pivot row\n", col, k)) ;
+		ASSERT (!NON_PIVOTAL_COL (col)) ;
+		Upos [col] = EMPTY ;
+		Lip [col] = EMPTY ;
+		Lilen [col] = 0 ;
+		k++ ;
+	    }
+	}
+	ASSERT (k == n_col) ;
+    }
+
+    if (npiv < n_inner)
+    {
+	/* finalize the diagonal of U */
+	DEBUGm3 (("Diag of U is zero, "ID" to "ID"\n", npiv, n_inner-1)) ;
+	for (k = npiv ; k < n_inner ; k++)
+	{
+	    CLEAR (D [k]) ;
+	}
+    }
+
+    /* save the pattern of the last row of U */
+    if (Numeric->ulen > 0)
+    {
+	DEBUGm3 (("Last row of U is not empty\n")) ;
+	Numeric->Upattern = Work->Upattern ;
+	Work->Upattern = (Int *) NULL ;
+    }
+
+    DEBUG2 (("Nnzpiv: "ID"  npiv "ID"\n", Numeric->nnzpiv, npiv)) ;
+    ASSERT (Numeric->nnzpiv <= npiv) ;
+    if (Numeric->nnzpiv < n_inner && !SCALAR_IS_NAN (Numeric->min_udiag))
+    {
+	/* the rest of the diagonal is zero, so min_udiag becomes 0,
+	 * unless it is already NaN. */
+	Numeric->min_udiag = 0.0 ;
+    }
+
+    /* ---------------------------------------------------------------------- */
+    /* size n_row, n_col workspaces that can be used here: */
+    /* ---------------------------------------------------------------------- */
+
+    Frpos = Work->Frpos ;	/* of size n_row+1 */
+    Fcpos = Work->Fcpos ;	/* of size n_col+1 */
+    Wp = Work->Wp ;		/* of size MAX(n_row,n_col)+1 */
+    /* Work->Upattern ;		cannot be used (in Numeric) */
+    Wr = Work->Lpattern ;	/* of size n_row+1 */
+    Wc = Work->Wrp ;		/* of size n_col+1 or bigger */
+
+    /* ---------------------------------------------------------------------- */
+    /* construct Rperm from inverse permutations */
+    /* ---------------------------------------------------------------------- */
+
+    /* use Frpos for temporary copy of inverse row permutation [ */
+
+    for (pivrow = 0 ; pivrow < n_row ; pivrow++)
+    {
+	k = Rperm [pivrow] ;
+	ASSERT (k < 0) ;
+	k = ONES_COMPLEMENT (k) ;
+	ASSERT (k >= 0 && k < n_row) ;
+	Wp [k] = pivrow ;
+	Frpos [pivrow] = k ;
+    }
+    for (k = 0 ; k < n_row ; k++)
+    {
+	Rperm [k] = Wp [k] ;
+    }
+
+    /* ---------------------------------------------------------------------- */
+    /* construct Cperm from inverse permutation */
+    /* ---------------------------------------------------------------------- */
+
+    /* use Fcpos for temporary copy of inverse column permutation [ */
+
+    for (pivcol = 0 ; pivcol < n_col ; pivcol++)
+    {
+	k = Cperm [pivcol] ;
+	ASSERT (k < 0) ;
+	k = ONES_COMPLEMENT (k) ;
+	ASSERT (k >= 0 && k < n_col) ;
+	Wp [k] = pivcol ;
+	/* save a copy of the inverse column permutation in Fcpos */
+	Fcpos [pivcol] = k ;
+    }
+    for (k = 0 ; k < n_col ; k++)
+    {
+	Cperm [k] = Wp [k] ;
+    }
+
+#ifndef NDEBUG
+    for (k = 0 ; k < n_col ; k++)
+    {
+	col = Cperm [k] ;
+	ASSERT (col >= 0 && col < n_col) ;
+	ASSERT (Fcpos [col] == k) ;		/* col is the kth pivot */
+    }
+    for (k = 0 ; k < n_row ; k++)
+    {
+	row = Rperm [k] ;
+	ASSERT (row >= 0 && row < n_row) ;
+	ASSERT (Frpos [row] == k) ;		/* row is the kth pivot */
+    }
+#endif
+
+#ifndef NDEBUG
+    UMF_dump_lu (Numeric) ;
+#endif
+
+    /* ---------------------------------------------------------------------- */
+    /* permute Lpos, Upos, Lilen, Lip, Uilen, and Uip */
+    /* ---------------------------------------------------------------------- */
+
+    for (k = 0 ; k < npiv ; k++)
+    {
+	pivrow = Rperm [k] ;
+	Wr [k] = Uilen [pivrow] ;
+	Wp [k] = Uip [pivrow] ;
+    }
+
+    for (k = 0 ; k < npiv ; k++)
+    {
+	Uilen [k] = Wr [k] ;
+	Uip [k] = Wp [k] ;
+    }
+
+    for (k = 0 ; k < npiv ; k++)
+    {
+	pivrow = Rperm [k] ;
+	Wp [k] = Lpos [pivrow] ;
+    }
+
+    for (k = 0 ; k < npiv ; k++)
+    {
+	Lpos [k] = Wp [k] ;
+    }
+
+    for (k = 0 ; k < npiv ; k++)
+    {
+	pivcol = Cperm [k] ;
+	Wc [k] = Lilen [pivcol] ;
+	Wp [k] = Lip [pivcol] ;
+    }
+
+    for (k = 0 ; k < npiv ; k++)
+    {
+	Lilen [k] = Wc [k] ;
+	Lip [k] = Wp [k] ;
+    }
+
+    for (k = 0 ; k < npiv ; k++)
+    {
+	pivcol = Cperm [k] ;
+	Wp [k] = Upos [pivcol] ;
+    }
+
+    for (k = 0 ; k < npiv ; k++)
+    {
+	Upos [k] = Wp [k] ;
+    }
+
+    /* ---------------------------------------------------------------------- */
+    /* terminate the last Uchain and last Lchain */
+    /* ---------------------------------------------------------------------- */
+
+    Upos [npiv] = EMPTY ;
+    Lpos [npiv] = EMPTY ;
+    Uip [npiv] = EMPTY ;
+    Lip [npiv] = EMPTY ;
+    Uilen [npiv] = 0 ;
+    Lilen [npiv] = 0 ;
+
+    /* ---------------------------------------------------------------------- */
+    /* convert U to the new pivot order */
+    /* ---------------------------------------------------------------------- */
+
+    n1 = Symbolic->n1 ;
+
+    for (k = 0 ; k < n1 ; k++)
+    {
+	/* this is a singleton row of U */
+	ulen = Uilen [k] ;
+	DEBUG4 (("K "ID" New U.  ulen "ID" Singleton 1\n", k, ulen)) ;
+	if (ulen > 0)
+	{
+	    up = Uip [k] ;
+	    ip = (Int *) (Numeric->Memory + up) ;
+	    for (i = 0 ; i < ulen ; i++)
+	    {
+		col = *ip ;
+		DEBUG4 ((" old col "ID" new col "ID"\n", col, Fcpos [col]));
+		ASSERT (col >= 0 && col < n_col) ;
+		*ip++ = Fcpos [col] ;
+	    }
+	}
+    }
+
+    for (k = n1 ; k < npiv ; k++)
+    {
+	up = Uip [k] ;
+	if (up < 0)
+	{
+	    /* this is the start of a new Uchain (with a pattern) */
+	    ulen = Uilen [k] ;
+	    DEBUG4 (("K "ID" New U.  ulen "ID" End_Uchain 1\n", k, ulen)) ;
+	    if (ulen > 0)
+	    {
+		up = -up ;
+		ip = (Int *) (Numeric->Memory + up) ;
+		for (i = 0 ; i < ulen ; i++)
+		{
+		    col = *ip ;
+		    DEBUG4 ((" old col "ID" new col "ID"\n", col, Fcpos [col]));
+		    ASSERT (col >= 0 && col < n_col) ;
+		    *ip++ = Fcpos [col] ;
+		}
+	    }
+	}
+    }
+
+    ulen = Numeric->ulen ;
+    if (ulen > 0)
+    {
+	/* convert last pivot row of U to the new pivot order */
+	DEBUG4 (("K "ID" (last)\n", k)) ;
+	for (i = 0 ; i < ulen ; i++)
+	{
+	    col = Numeric->Upattern [i] ;
+	    DEBUG4 (("    old col "ID" new col "ID"\n", col, Fcpos [col])) ;
+	    Numeric->Upattern [i] = Fcpos [col] ;
+	}
+    }
+
+    /* Fcpos no longer needed ] */
+
+    /* ---------------------------------------------------------------------- */
+    /* convert L to the new pivot order */
+    /* ---------------------------------------------------------------------- */
+
+    for (k = 0 ; k < n1 ; k++)
+    {
+	llen = Lilen [k] ;
+	DEBUG4 (("K "ID" New L.  llen "ID" Singleton col\n", k, llen)) ;
+	if (llen > 0)
+	{
+	    lp = Lip [k] ;
+	    ip = (Int *) (Numeric->Memory + lp) ;
+	    for (i = 0 ; i < llen ; i++)
+	    {
+		row = *ip ;
+		DEBUG4 (("    old row "ID" new row "ID"\n", row, Frpos [row])) ;
+		ASSERT (row >= 0 && row < n_row) ;
+		*ip++ = Frpos [row] ;
+	    }
+	}
+    }
+
+    for (k = n1 ; k < npiv ; k++)
+    {
+	llen = Lilen [k] ;
+	DEBUG4 (("K "ID" New L.  llen "ID" \n", k, llen)) ;
+	if (llen > 0)
+	{
+	    lp = Lip [k] ;
+	    if (lp < 0)
+	    {
+		/* this starts a new Lchain */
+		lp = -lp ;
+	    }
+	    ip = (Int *) (Numeric->Memory + lp) ;
+	    for (i = 0 ; i < llen ; i++)
+	    {
+		row = *ip ;
+		DEBUG4 (("    old row "ID" new row "ID"\n", row, Frpos [row])) ;
+		ASSERT (row >= 0 && row < n_row) ;
+		*ip++ = Frpos [row] ;
+	    }
+	}
+    }
+
+    /* Frpos no longer needed ] */
+
+    /* ---------------------------------------------------------------------- */
+    /* combine symbolic and numeric permutations */
+    /* ---------------------------------------------------------------------- */
+
+    Cperm_init = Symbolic->Cperm_init ;
+    Rperm_init = Symbolic->Rperm_init ;
+
+    for (k = 0 ; k < n_row ; k++)
+    {
+	Rperm [k] = Rperm_init [Rperm [k]] ;
+    }
+
+    for (k = 0 ; k < n_col ; k++)
+    {
+	Cperm [k] = Cperm_init [Cperm [k]] ;
+    }
+
+    /* Work object will be freed immediately upon return (to UMF_kernel */
+    /* and then to UMFPACK_numeric). */
+}