/* ========================================================================== */
/* === UMFPACK_report_numeric =============================================== */
/* ========================================================================== */

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

/*
    User-callable.  Prints the Numeric object.
    See umfpack_report_numeric.h for details.

    Dynamic memory usage:  Allocates a size n*sizeof(Int) workspace via a single
    call to UMF_malloc and then frees all of it via UMF_free on return.  The
    workspace is not allocated if an early error return occurs before the
    workspace is needed.
*/

#include "umf_internal.h"
#include "umf_valid_numeric.h"
#include "umf_report_perm.h"
#include "umf_report_vector.h"
#include "umf_malloc.h"
#include "umf_free.h"


PRIVATE Int report_L
(
    NumericType *Numeric,
    Int Pattern [ ],
    Int prl
) ;


PRIVATE Int report_U
(
    NumericType *Numeric,
    Int Pattern [ ],
    Int prl
) ;

/* ========================================================================== */
/* === UMFPACK_report_numeric =============================================== */
/* ========================================================================== */

GLOBAL Int UMFPACK_report_numeric
(
    void *NumericHandle,
    const double Control [UMFPACK_CONTROL]
)
{
    Int prl, *W, nn, n_row, n_col, n_inner, num_fixed_size, numeric_size,
	npiv ;
    NumericType *Numeric ;

    prl = GET_CONTROL (UMFPACK_PRL, UMFPACK_DEFAULT_PRL) ;

    if (prl <= 2)
    {
	return (UMFPACK_OK) ;
    }

    PRINTF (("Numeric object:  ")) ;

    Numeric = (NumericType *) NumericHandle ;
    if (!UMF_valid_numeric (Numeric))
    {
	PRINTF (("ERROR: LU factors invalid\n\n")) ;
	return (UMFPACK_ERROR_invalid_Numeric_object) ;
    }

    n_row = Numeric->n_row ;
    n_col = Numeric->n_col ;
    nn = MAX (n_row, n_col) ;
    n_inner = MIN (n_row, n_col) ;
    npiv = Numeric->npiv ;

    DEBUG1 (("n_row "ID" n_col "ID" nn "ID" n_inner "ID" npiv "ID"\n",
	n_row, n_col, nn, n_inner, npiv)) ;

    /* size of Numeric object, except Numeric->Memory and Numeric->Upattern */
    /* see also UMF_set_stats */
    num_fixed_size =
	UNITS (NumericType, 1)		/* Numeric structure */
	+ UNITS (Entry, n_inner+1)	/* D */
	+ UNITS (Int, n_row+1)		/* Rperm */
	+ UNITS (Int, n_col+1)		/* Cperm */
	+ 6 * UNITS (Int, npiv+1)	/* Lpos, Uilen, Uip, Upos, Lilen, Lip */
	+ ((Numeric->scale != UMFPACK_SCALE_NONE) ?
		UNITS (Entry, n_row) : 0) ; /* Rs */

    DEBUG1 (("num fixed size: "ID"\n", num_fixed_size)) ;
    DEBUG1 (("Numeric->size "ID"\n", Numeric->size)) ;
    DEBUG1 (("ulen units "ID"\n", UNITS (Int, Numeric->ulen))) ;

    /* size of Numeric->Memory is Numeric->size */
    /* size of Numeric->Upattern is Numeric->ulen */
    numeric_size = num_fixed_size + Numeric->size
	+ UNITS (Int, Numeric->ulen) ;

    DEBUG1 (("numeric total size "ID"\n", numeric_size)) ;

    if (prl >= 4)
    {
	PRINTF (("\n    n_row: "ID"  n_col: "ID"\n", n_row, n_col)) ;

	PRINTF (("    relative pivot tolerance used:              %g\n",
	    Numeric->relpt)) ;
	PRINTF (("    relative symmetric pivot tolerance used:    %g\n",
	    Numeric->relpt2)) ;

	PRINTF (("    matrix scaled: ")) ;
	if (Numeric->scale == UMFPACK_SCALE_NONE)
	{
	    PRINTF (("no")) ;
	}
	else if (Numeric->scale == UMFPACK_SCALE_SUM)
	{
	    PRINTF (("yes (divided each row by sum abs value in each row)\n")) ;
	    PRINTF (("    minimum sum (abs (rows of A)):              %.5e\n",
		Numeric->rsmin)) ;
	    PRINTF (("    maximum sum (abs (rows of A)):              %.5e",
		Numeric->rsmax)) ;
	}
	else if (Numeric->scale == UMFPACK_SCALE_MAX)
	{
	    PRINTF (("yes (divided each row by max abs value in each row)\n")) ;
	    PRINTF (("    minimum max (abs (rows of A)):              %.5e\n",
		Numeric->rsmin)) ;
	    PRINTF (("    maximum max (abs (rows of A)):              %.5e",
		Numeric->rsmax)) ;
	}
	PRINTF (("\n")) ;

	PRINTF (("    initial allocation parameter used:          %g\n",
	    Numeric->alloc_init)) ;
	PRINTF (("    frontal matrix allocation parameter used:   %g\n",
	    Numeric->front_alloc_init)) ;
	PRINTF (("    final total size of Numeric object (Units): "ID"\n",
	    numeric_size)) ;
	PRINTF (("    final total size of Numeric object (MBytes): %.1f\n",
	    MBYTES (numeric_size))) ;
	PRINTF (("    peak size of variable-size part (Units):    "ID"\n",
	    Numeric->max_usage)) ;
	PRINTF (("    peak size of variable-size part (MBytes):   %.1f\n",
	    MBYTES (Numeric->max_usage))) ;
	PRINTF (("    largest actual frontal matrix size:         "ID"\n",
	    Numeric->maxfrsize)) ;
	PRINTF (("    memory defragmentations:                    "ID"\n",
	    Numeric->ngarbage)) ;
	PRINTF (("    memory reallocations:                       "ID"\n",
	    Numeric->nrealloc)) ;
	PRINTF (("    costly memory reallocations:                "ID"\n",
	    Numeric->ncostly)) ;
	PRINTF (("    entries in compressed pattern (L and U):    "ID"\n",
	    Numeric->isize)) ;
	PRINTF (("    number of nonzeros in L (excl diag):        "ID"\n",
	    Numeric->lnz)) ;
	PRINTF (("    number of entries stored in L (excl diag):  "ID"\n",
	    Numeric->nLentries)) ;
	PRINTF (("    number of nonzeros in U (excl diag):        "ID"\n",
	    Numeric->unz)) ;
	PRINTF (("    number of entries stored in U (excl diag):  "ID"\n",
	    Numeric->nUentries)) ;
	PRINTF (("    factorization floating-point operations:    %g\n",
	    Numeric->flops)) ;
	PRINTF (("    number of nonzeros on diagonal of U:        "ID"\n",
	    Numeric->nnzpiv)) ;
	PRINTF (("    min abs. value on diagonal of U:            %.5e\n",
	    Numeric->min_udiag)) ;
	PRINTF (("    max abs. value on diagonal of U:            %.5e\n",
	    Numeric->max_udiag)) ;
	PRINTF (("    reciprocal condition number estimate:       %.2e\n",
	    Numeric->rcond)) ;
    }

    W = (Int *) UMF_malloc (nn, sizeof (Int)) ;
    if (!W)
    {
	PRINTF ((" ERROR: out of memory to check Numeric object\n\n")) ;
	return (UMFPACK_ERROR_out_of_memory) ;
    }

    if (Numeric->Rs)
    {
#ifndef NRECIPROCAL
	if (Numeric->do_recip)
	{
	    PRINTF4 (("\nScale factors applied via multiplication\n")) ;
	}
	else
#endif
	{
	    PRINTF4 (("\nScale factors applied via division\n")) ;
	}
	PRINTF4 (("Scale factors, Rs: ")) ;
	(void) UMF_report_vector (n_row, Numeric->Rs, (double *) NULL,
	    prl, FALSE, TRUE) ;
    }
    else
    {
	PRINTF4 (("Scale factors, Rs: (not present)\n")) ;
    }

    PRINTF4 (("\nP: row ")) ;
    if (UMF_report_perm (n_row, Numeric->Rperm, W, prl, 0) != UMFPACK_OK)
    {
	(void) UMF_free ((void *) W) ;
	return (UMFPACK_ERROR_invalid_Numeric_object) ;
    }

    PRINTF4 (("\nQ: column ")) ;
    if (UMF_report_perm (n_col, Numeric->Cperm, W, prl, 0) != UMFPACK_OK)
    {
	(void) UMF_free ((void *) W) ;
	return (UMFPACK_ERROR_invalid_Numeric_object) ;
    }

    if (!report_L (Numeric, W, prl))
    {
	(void) UMF_free ((void *) W) ;
	PRINTF ((" ERROR: L factor invalid\n\n")) ;
	return (UMFPACK_ERROR_invalid_Numeric_object) ;
    }

    if (!report_U (Numeric, W, prl))
    {
	(void) UMF_free ((void *) W) ;
	PRINTF ((" ERROR: U factor invalid\n\n")) ;
	return (UMFPACK_ERROR_invalid_Numeric_object) ;
    }

    /* The diagonal of U is in "merged" (Entry) form, not "split" form. */
    PRINTF4 (("\ndiagonal of U: ")) ;
    (void) UMF_report_vector (n_inner, (double *) Numeric->D, (double *) NULL,
	prl, FALSE, FALSE) ;

    (void) UMF_free ((void *) W) ;

    PRINTF4 (("    Numeric object:  ")) ;
    PRINTF (("OK\n\n")) ;
    return (UMFPACK_OK) ;
}


/* ========================================================================== */
/* === report_L ============================================================= */
/* ========================================================================== */

PRIVATE Int report_L
(
    NumericType *Numeric,
    Int Pattern [ ],
    Int prl
)
{
    Int k, deg, *ip, j, row, n_row, *Lpos, *Lilen, valid, k1,
	*Lip, newLchain, llen, prl1, pos, lp, p, npiv, n1, *Li ;
    Entry *xp, *Lval ;

    /* ---------------------------------------------------------------------- */

    ASSERT (prl >= 3) ;

    n_row = Numeric->n_row ;
    npiv = Numeric->npiv ;
    n1 = Numeric->n1 ;
    Lpos = Numeric->Lpos ;
    Lilen = Numeric->Lilen ;
    Lip = Numeric->Lip ;
    prl1 = prl ;
    deg = 0 ;

    PRINTF4 ((
    "\nL in Numeric object, in column-oriented compressed-pattern form:\n"
    "    Diagonal entries are all equal to 1.0 (not stored)\n")) ;

    ASSERT (Pattern != (Int *) NULL) ;

    /* ---------------------------------------------------------------------- */
    /* print L */
    /* ---------------------------------------------------------------------- */

    k1 = 12 ;

    /* ---------------------------------------------------------------------- */
    /* print the singleton columns of L */
    /* ---------------------------------------------------------------------- */

    for (k = 0 ; k < n1 ; k++)
    {
	if (k1 > 0)
	{
	    prl = prl1 ;
	}
	lp = Lip [k] ;
	deg = Lilen [k] ;
	Li = (Int *) (Numeric->Memory + lp) ;
	lp += UNITS (Int, deg) ;
	Lval = (Entry *) (Numeric->Memory + lp) ;
	if (k1-- > 0)
	{
	    prl = prl1 ;
	}
	else if (prl == 4)
	{
	    PRINTF (("    ...\n")) ;
	    prl-- ;
	}
	PRINTF4 (("\n    column "ID":", INDEX (k))) ;
	PRINTF4 (("  length "ID".\n", deg)) ;
	for (j = 0 ; j < deg ; j++)
	{
	    row = Li [j] ;
	    PRINTF4 (("\trow "ID" : ", INDEX (row))) ;
	    if (prl >= 4) PRINT_ENTRY (Lval [j]) ;
	    if (row <= k || row >= n_row)
	    {
		return (FALSE) ;
	    }
	    PRINTF4 (("\n")) ;
	    /* truncate printout, but continue to check L */
	    if (prl == 4 && j == 9 && deg > 10)
	    {
		PRINTF (("\t...\n")) ;
		prl-- ;
	    }
	}
    }

    /* ---------------------------------------------------------------------- */
    /* print the regular columns of L */
    /* ---------------------------------------------------------------------- */

    for (k = n1 ; k < npiv ; k++)
    {
	/* if prl is 4, print the first 10 entries of the first 10 columns */
	if (k1 > 0)
	{
	    prl = prl1 ;
	}

	lp = Lip [k] ;
	newLchain = (lp < 0) ;
	if (newLchain)
	{
	    lp = -lp ;
	    deg = 0 ;
	}

	if (k1-- > 0)
	{
	    prl = prl1 ;
	}
	else if (prl == 4)
	{
	    PRINTF (("    ...\n")) ;
	    prl-- ;
	}

	PRINTF4 (("\n    column "ID":", INDEX (k))) ;

	/* ------------------------------------------------------------------ */
	/* make column of L in Pattern [0..deg-1] */
	/* ------------------------------------------------------------------ */

	/* remove pivot row */
	pos = Lpos [k] ;
	if (pos != EMPTY)
	{
	    PRINTF4 (("  remove row "ID" at position "ID".",
		INDEX (Pattern [pos]), INDEX (pos))) ;
	    valid = (!newLchain) && (deg > 0) && (pos < deg) && (pos >= 0)
		&& (Pattern [pos] == k) ;
	    if (!valid)
	    {
		return (FALSE) ;
	    }
	    Pattern [pos] = Pattern [--deg] ;
	}

	/* concatenate the pattern */
	llen = Lilen [k] ;
	if (llen < 0)
	{
	    return (FALSE) ;
	}
	p = lp + UNITS (Int, llen) ;
	xp = (Entry *) (Numeric->Memory + p) ;
	if ((llen > 0 || deg > 0)
	    && (p + (Int) UNITS (Entry, deg) > Numeric->size))
	{
	    return (FALSE) ;
	}
	if (llen > 0)
	{
	    PRINTF4 (("  add "ID" entries.", llen)) ;
	    ip = (Int *) (Numeric->Memory + lp) ;
	    for (j = 0 ; j < llen ; j++)
	    {
		Pattern [deg++] = *ip++ ;
	    }
	}

	/* ------------------------------------------------------------------ */
	/* print column k of L */
	/* ------------------------------------------------------------------ */

	PRINTF4 (("  length "ID".", deg)) ;
	if (newLchain)
	{
	    PRINTF4 (("  Start of Lchain.")) ;
	}
	PRINTF4 (("\n")) ;

	for (j = 0 ; j < deg ; j++)
	{
	    row = Pattern [j] ;
	    PRINTF4 (("\trow "ID" : ", INDEX (row))) ;
	    if (prl >= 4) PRINT_ENTRY (*xp) ;
	    if (row <= k || row >= n_row)
	    {
		return (FALSE) ;
	    }
	    PRINTF4 (("\n")) ;
	    xp++ ;
	    /* truncate printout, but continue to check L */
	    if (prl == 4 && j == 9 && deg > 10)
	    {
		PRINTF (("\t...\n")) ;
		prl-- ;
	    }
	}
    }

    PRINTF4 (("\n")) ;
    return (TRUE) ;
}


/* ========================================================================== */
/* === report_U ============================================================= */
/* ========================================================================== */

PRIVATE Int report_U
(
    NumericType *Numeric,
    Int Pattern [ ],
    Int prl
)
{
    /* ---------------------------------------------------------------------- */

    Int k, deg, j, *ip, col, *Upos, *Uilen, k1, prl1, pos,
	*Uip, n_col, ulen, p, newUchain, up, npiv, n1, *Ui ;
    Entry *xp, *Uval ;

    /* ---------------------------------------------------------------------- */

    ASSERT (prl >= 3) ;

    n_col = Numeric->n_col ;
    npiv = Numeric->npiv ;
    n1 = Numeric->n1 ;
    Upos = Numeric->Upos ;
    Uilen = Numeric->Uilen ;
    Uip = Numeric->Uip ;
    prl1 = prl ;

    PRINTF4 ((
    "\nU in Numeric object, in row-oriented compressed-pattern form:\n"
    "    Diagonal is stored separately.\n")) ;

    ASSERT (Pattern != (Int *) NULL) ;

    k1 = 12 ;

    /* ---------------------------------------------------------------------- */
    /* print the sparse part of U */
    /* ---------------------------------------------------------------------- */

    deg = Numeric->ulen ;
    if (deg > 0)
    {
	/* make last pivot row of U (singular matrices only) */
	for (j = 0 ; j < deg ; j++)
	{
	    Pattern [j] = Numeric->Upattern [j] ;
	}
    }

    PRINTF4 (("\n    row "ID":  length "ID".  End of Uchain.\n", INDEX (npiv-1),
	deg)) ;

    for (k = npiv-1 ; k >= n1 ; k--)
    {

	/* ------------------------------------------------------------------ */
	/* print row k of U */
	/* ------------------------------------------------------------------ */

	/* if prl is 3, print the first 10 entries of the first 10 columns */
	if (k1 > 0)
	{
	    prl = prl1 ;
	}

	up = Uip [k] ;
	ulen = Uilen [k] ;
	if (ulen < 0)
	{
	    return (FALSE) ;
	}
	newUchain = (up < 0) ;
	if (newUchain)
	{
	    up = -up ;
	    p = up + UNITS (Int, ulen) ;
	}
	else
	{
	    p = up ;
	}
	xp = (Entry *) (Numeric->Memory + p) ;
	if (deg > 0 && (p + (Int) UNITS (Entry, deg) > Numeric->size))
	{
	    return (FALSE) ;
	}
	for (j = 0 ; j < deg ; j++)
	{
	    col = Pattern [j] ;
	    PRINTF4 (("\tcol "ID" :", INDEX (col))) ;
	    if (prl >= 4) PRINT_ENTRY (*xp) ;
	    if (col <= k || col >= n_col)
	    {
		return (FALSE) ;
	    }
	    PRINTF4 (("\n")) ;
	    xp++ ;
	    /* truncate printout, but continue to check U */
	    if (prl == 4 && j == 9 && deg > 10)
	    {
		PRINTF (("\t...\n")) ;
		prl-- ;
	    }
	}

	/* ------------------------------------------------------------------ */
	/* make row k-1 of U in Pattern [0..deg-1] */
	/* ------------------------------------------------------------------ */

	if (k1-- > 0)
	{
	    prl = prl1 ;
	}
	else if (prl == 4)
	{
	    PRINTF (("    ...\n")) ;
	    prl-- ;
	}

	if (k > 0)
	{
	    PRINTF4 (("\n    row "ID":  ", INDEX (k-1))) ;
	}

	if (newUchain)
	{
	    /* next row is a new Uchain */
	    if (k > 0)
	    {
		deg = ulen ;
		PRINTF4 (("length "ID".  End of Uchain.\n", deg)) ;
		if (up + (Int) UNITS (Int, ulen) > Numeric->size)
		{
		    return (FALSE) ;
		}
		ip = (Int *) (Numeric->Memory + up) ;
		for (j = 0 ; j < deg ; j++)
		{
		    Pattern [j] = *ip++ ;
		}
	    }
	}
	else
	{
	    if (ulen > 0)
	    {
		PRINTF4 (("remove "ID" entries.  ", ulen)) ;
	    }
	    deg -= ulen ;
	    if (deg < 0)
	    {
		return (FALSE) ;
	    }
	    pos = Upos [k] ;
	    if (pos != EMPTY)
	    {
		/* add the pivot column */
		PRINTF4 (("add column "ID" at position "ID".  ",
		    INDEX (k), INDEX (pos))) ;
		if (pos < 0 || pos > deg)
		{
		    return (FALSE) ;
		}
		Pattern [deg++] = Pattern [pos] ;
		Pattern [pos] = k ;
	    }
	    PRINTF4 (("length "ID".\n", deg)) ;
	}
    }

    /* ---------------------------------------------------------------------- */
    /* print the singleton rows of U */
    /* ---------------------------------------------------------------------- */

    for (k = n1 - 1 ; k >= 0 ; k--)
    {
	if (k1 > 0)
	{
	    prl = prl1 ;
	}
	up = Uip [k] ;
	deg = Uilen [k] ;
	Ui = (Int *) (Numeric->Memory + up) ;
	up += UNITS (Int, deg) ;
	Uval = (Entry *) (Numeric->Memory + up) ;
	if (k1-- > 0)
	{
	    prl = prl1 ;
	}
	else if (prl == 4)
	{
	    PRINTF (("    ...\n")) ;
	    prl-- ;
	}
	PRINTF4 (("\n    row "ID":", INDEX (k))) ;
	PRINTF4 (("  length "ID".\n", deg)) ;
	for (j = 0 ; j < deg ; j++)
	{
	    col = Ui [j] ;
	    PRINTF4 (("\tcol "ID" : ", INDEX (col))) ;
	    if (prl >= 4) PRINT_ENTRY (Uval [j]) ;
	    if (col <= k || col >= n_col)
	    {
		return (FALSE) ;
	    }
	    PRINTF4 (("\n")) ;
	    /* truncate printout, but continue to check U */
	    if (prl == 4 && j == 9 && deg > 10)
	    {
		PRINTF (("\t...\n")) ;
		prl-- ;
	    }
	}
    }

    prl = prl1 ;
    PRINTF4 (("\n")) ;
    return (TRUE) ;
}