/* ========================================================================== */
/* === UMF_kernel =========================================================== */
/* ========================================================================== */

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

/*
    Primary factorization routine.   Called by UMFPACK_numeric.
    Returns:
	UMFPACK_OK if successful,
	UMFPACK_ERROR_out_of_memory if out of memory, or
	UMFPACK_ERROR_different_pattern if pattern of matrix (Ap and/or Ai)
	   has changed since the call to UMFPACK_*symbolic.
*/

#include "umf_internal.h"
#include "umf_kernel_init.h"
#include "umf_init_front.h"
#include "umf_start_front.h"
#include "umf_assemble.h"
#include "umf_scale_column.h"
#include "umf_local_search.h"
#include "umf_create_element.h"
#include "umf_extend_front.h"
#include "umf_blas3_update.h"
#include "umf_store_lu.h"
#include "umf_kernel_wrapup.h"

/* perform an action, and return if out of memory */
#define DO(action) { if (! (action)) { return (UMFPACK_ERROR_out_of_memory) ; }}

GLOBAL Int UMF_kernel
(
    const Int Ap [ ],
    const Int Ai [ ],
    const double Ax [ ],
#ifdef COMPLEX
    const double Az [ ],
#endif
    NumericType *Numeric,
    WorkType *Work,
    SymbolicType *Symbolic
)
{

    /* ---------------------------------------------------------------------- */
    /* local variables */
    /* ---------------------------------------------------------------------- */

    Int j, f1, f2, chain, nchains, *Chain_start, status, fixQ, evaporate,
	*Front_npivcol, jmax, nb, drop ;

    /* ---------------------------------------------------------------------- */
    /* initialize memory space and load the matrix. Optionally scale. */
    /* ---------------------------------------------------------------------- */

    if (!UMF_kernel_init (Ap, Ai, Ax,
#ifdef COMPLEX
	Az,
#endif
	Numeric, Work, Symbolic))
    {
	/* UMF_kernel_init is guaranteed to succeed, since UMFPACK_numeric */
	/* either allocates enough space or if not, UMF_kernel does not get */
	/* called.  So running out of memory here is a fatal error, and means */
	/* that the user changed Ap and/or Ai since the call to */
	/* UMFPACK_*symbolic. */
	DEBUGm4 (("kernel init failed\n")) ;
	return (UMFPACK_ERROR_different_pattern) ;
    }

    /* ---------------------------------------------------------------------- */
    /* get the symbolic factorization */
    /* ---------------------------------------------------------------------- */

    nchains = Symbolic->nchains ;
    Chain_start = Symbolic->Chain_start ;
    Front_npivcol = Symbolic->Front_npivcol ;
    nb = Symbolic->nb ;
    fixQ = Symbolic->fixQ ;
    drop = Numeric->droptol > 0.0 ;

#ifndef NDEBUG
    for (chain = 0 ; chain < nchains ; chain++)
    {
	Int i ;
	f1 = Chain_start [chain] ;
	f2 = Chain_start [chain+1] - 1 ;
	DEBUG1 (("\nCHain: "ID" start "ID" end "ID"\n", chain, f1, f2)) ;
	for (i = f1 ; i <= f2 ; i++)
	{
	    DEBUG1 (("Front "ID", npivcol "ID"\n", i, Front_npivcol [i])) ;
	}
    }
#endif

    /* ---------------------------------------------------------------------- */
    /* factorize each chain of frontal matrices */
    /* ---------------------------------------------------------------------- */

    for (chain = 0 ; chain < nchains ; chain++)
    {
	f1 = Chain_start [chain] ;
	f2 = Chain_start [chain+1] - 1 ;

	/* ------------------------------------------------------------------ */
	/* get the initial frontal matrix size for this chain */
	/* ------------------------------------------------------------------ */

	DO (UMF_start_front (chain, Numeric, Work, Symbolic)) ;

	/* ------------------------------------------------------------------ */
	/* factorize each front in the chain */
	/* ------------------------------------------------------------------ */

	for (Work->frontid = f1 ; Work->frontid <= f2 ; Work->frontid++)
	{

	    /* -------------------------------------------------------------- */
	    /* Initialize the pivot column candidate set  */
	    /* -------------------------------------------------------------- */

	    Work->ncand = Front_npivcol [Work->frontid] ;
	    Work->lo = Work->nextcand ;
	    Work->hi = Work->nextcand + Work->ncand - 1 ;
	    jmax = MIN (MAX_CANDIDATES, Work->ncand) ;
	    DEBUGm1 ((">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> Starting front "
		ID", npivcol: "ID"\n", Work->frontid, Work->ncand)) ;
	    if (fixQ)
	    {
		/* do not modify the column order */
		jmax = 1 ;
	    }
	    DEBUGm1 (("Initial candidates: ")) ;
	    for (j = 0 ; j < jmax ; j++)
	    {
		DEBUGm1 ((" "ID, Work->nextcand)) ;
		ASSERT (Work->nextcand <= Work->hi) ;
		Work->Candidates [j] = Work->nextcand++ ;
	    }
	    Work->nCandidates = jmax ;
	    DEBUGm1 (("\n")) ;

	    /* -------------------------------------------------------------- */
	    /* Assemble and factorize the current frontal matrix */
	    /* -------------------------------------------------------------- */

	    while (Work->ncand > 0)
	    {

		/* ---------------------------------------------------------- */
		/* get the pivot row and column */
		/* ---------------------------------------------------------- */

		status = UMF_local_search (Numeric, Work, Symbolic) ;
		if (status == UMFPACK_ERROR_different_pattern)
		{
		    /* :: pattern change detected in umf_local_search :: */
		    /* input matrix has changed since umfpack_*symbolic */
		    DEBUGm4 (("local search failed\n")) ;
		    return (UMFPACK_ERROR_different_pattern) ;
		}
		if (status == UMFPACK_WARNING_singular_matrix)
		{
		    /* no pivot found, discard and try again */
		    continue ;
		}

		/* ---------------------------------------------------------- */
		/* update if front not extended or too many zeros in L,U */
		/* ---------------------------------------------------------- */

		if (Work->do_update)
		{
		    UMF_blas3_update (Work) ;
		    if (drop)
		    {
			DO (UMF_store_lu_drop (Numeric, Work)) ;
		    }
		    else
		    {
			DO (UMF_store_lu (Numeric, Work)) ;
		    }
		}

		/* ---------------------------------------------------------- */
		/* extend the frontal matrix, or start a new one */
		/* ---------------------------------------------------------- */

		if (Work->do_extend)
		{
		    /* extend the current front */
		    DO (UMF_extend_front (Numeric, Work)) ;
		}
		else
		{
		    /* finish the current front (if any) and start a new one */
		    DO (UMF_create_element (Numeric, Work, Symbolic)) ;
		    DO (UMF_init_front (Numeric, Work)) ;
		}

		/* ---------------------------------------------------------- */
		/* Numerical & symbolic assembly into current frontal matrix */
		/* ---------------------------------------------------------- */

		if (fixQ)
		{
		    UMF_assemble_fixq (Numeric, Work) ;
		}
		else
		{
		    UMF_assemble (Numeric, Work) ;
		}

		/* ---------------------------------------------------------- */
		/* scale the pivot column */
		/* ---------------------------------------------------------- */

		UMF_scale_column (Numeric, Work) ;

		/* ---------------------------------------------------------- */
		/* Numerical update if enough pivots accumulated */
		/* ---------------------------------------------------------- */

		evaporate = Work->fnrows == 0 || Work->fncols == 0 ;
		if (Work->fnpiv >= nb || evaporate)
		{
		    UMF_blas3_update (Work) ;
		    if (drop)
		    {
			DO (UMF_store_lu_drop (Numeric, Work)) ;
		    }
		    else
		    {
			DO (UMF_store_lu (Numeric, Work)) ;
		    }

		}

		Work->pivrow_in_front = FALSE ;
		Work->pivcol_in_front = FALSE ;

		/* ---------------------------------------------------------- */
		/* If front is empty, evaporate it */
		/* ---------------------------------------------------------- */

		if (evaporate)
		{
		    /* This does not create an element, just evaporates it.
		     * It ensures that a front is not 0-by-c or r-by-0.  No
		     * memory is allocated, so it is guaranteed to succeed. */
		    (void) UMF_create_element (Numeric, Work, Symbolic) ;
		    Work->fnrows = 0 ;
		    Work->fncols = 0 ;
		}
	    }
	}

	/* ------------------------------------------------------------------
	 * Wrapup the current frontal matrix.  This is the last in a chain
	 * in the column elimination tree.  The next frontal matrix
	 * cannot overlap with the current one, which will be its sibling
	 * in the column etree.
	 * ------------------------------------------------------------------ */

	UMF_blas3_update (Work) ;
	if (drop)
	{
	    DO (UMF_store_lu_drop (Numeric, Work)) ;
	}
	else
	{
	    DO (UMF_store_lu (Numeric, Work)) ;
	}
	Work->fnrows_new = Work->fnrows ;
	Work->fncols_new = Work->fncols ;
	DO (UMF_create_element (Numeric, Work, Symbolic)) ;

	/* ------------------------------------------------------------------ */
	/* current front is now empty */
	/* ------------------------------------------------------------------ */

	Work->fnrows = 0 ;
	Work->fncols = 0 ;
    }

    /* ---------------------------------------------------------------------- */
    /* end the last Lchain and Uchain and finalize the LU factors */
    /* ---------------------------------------------------------------------- */

    UMF_kernel_wrapup (Numeric, Symbolic, Work) ;

    /* note that the matrix may be singular (this is OK) */
    return (UMFPACK_OK) ;
}