//------------------------------------------------------------------------------
// CHOLMOD/Supernodal/cholmod_super_numeric: supernodal Cholesky factorization
//------------------------------------------------------------------------------

// CHOLMOD/Supernodal Module.  Copyright (C) 2005-2023, Timothy A. Davis.
// All Rights Reserved.
// SPDX-License-Identifier: GPL-2.0+

//------------------------------------------------------------------------------

// Computes the Cholesky factorization of A+beta*I or A*F+beta*I.  Only the
// the lower triangular part of A+beta*I or A*F+beta*I is accessed.  The
// matrices A and F must already be permuted according to the fill-reduction
// permutation L->Perm.  cholmod_factorize is an "easy" wrapper for this code
// which applies that permutation.  beta is real.
//
// Symmetric case: A is a symmetric (lower) matrix.  F is not accessed.
// With a fill-reducing permutation, A(p,p) should be passed instead, where is
// p is L->Perm.
//
// Unsymmetric case: A is unsymmetric, and F must be present.  Normally, F=A'.
// With a fill-reducing permutation, A(p,f) and A(p,f)' should be passed as A
// and F, respectively, where f is a list of the subset of the columns of A.
//
// The input factorization L must be supernodal (L->is_super is TRUE).  It can
// either be symbolic or numeric.  In the first case, L has been analyzed by
// cholmod_analyze or cholmod_super_symbolic, but the matrix has not yet been
// numerically factorized.  The numerical values are allocated here and the
// factorization is computed.  In the second case, a prior matrix has been
// analyzed and numerically factorized, and a new matrix is being factorized.
// The numerical values of L are replaced with the new numerical factorization.
//
// L->is_ll is ignored, and set to TRUE.  This routine always computes an LL'
// factorization.  Supernodal LDL' factorization is not (yet) supported.
// FUTURE WORK: perform a supernodal LDL' factorization if L->is_ll is FALSE.
//
// If the matrix is not positive definite the routine returns TRUE, but sets
// Common->status to CHOLMOD_NOT_POSDEF and L->minor is set to the column at
// which the failure occurred.  The supernode containing the non-positive
// diagonal entry is set to zero (this includes columns to the left of L->minor
// in the same supernode), as are all subsequent supernodes.
//
// workspace: Flag (nrow), Head (nrow+1), Iwork (2*nrow + 5*nsuper).
//      Allocates temporary space of size L->maxcsize * sizeof(double)
//      (twice that for the complex/zomplex case).
//
// If L is supernodal symbolic on input, it is converted to a supernodal numeric
// factor on output, with an xtype of real if A is real, or complex if A is
// complex or zomplex.  If L is supernodal numeric on input, its xtype must
// match A (except that L can be complex and A zomplex).  The xtype of A and F
// must match.  The dtypes of all matrices must match.

#include "cholmod_internal.h"

#ifndef NGPL
#ifndef NSUPERNODAL

//------------------------------------------------------------------------------
// GPU templates: double and double complex cases only
//------------------------------------------------------------------------------

#define DOUBLE
#if defined ( CHOLMOD_HAS_CUDA ) && defined ( CHOLMOD_INT64 )
    #include "cholmod_gpu_kernels.h"
    #define REAL
    #include "../GPU/t_cholmod_gpu.c"
    #define COMPLEX
    #include "../GPU/t_cholmod_gpu.c"
#endif

//------------------------------------------------------------------------------
// CPU templates:  all cases
//------------------------------------------------------------------------------

#define DOUBLE
#define REAL
#include "t_cholmod_super_numeric_worker.c"
#define COMPLEX
#include "t_cholmod_super_numeric_worker.c"
#define ZOMPLEX
#include "t_cholmod_super_numeric_worker.c"

#undef  DOUBLE
#define SINGLE
#define REAL
#include "t_cholmod_super_numeric_worker.c"
#define COMPLEX
#include "t_cholmod_super_numeric_worker.c"
#define ZOMPLEX
#include "t_cholmod_super_numeric_worker.c"

//------------------------------------------------------------------------------
// cholmod_super_numeric
//------------------------------------------------------------------------------

// Returns TRUE if successful, or if the matrix is not positive definite.
// Returns FALSE if out of memory, inputs are invalid, or other fatal error
// occurs.

int CHOLMOD(super_numeric)
(
    // input:
    cholmod_sparse *A,  // matrix to factorize
    cholmod_sparse *F,  // F = A' or A(:,f)'
    double beta [2],    // beta*I is added to diagonal of matrix to factorize
    // input/output:
    cholmod_factor *L,  // factorization
    cholmod_common *Common
)
{

    //--------------------------------------------------------------------------
    // check inputs
    //--------------------------------------------------------------------------

    cholmod_dense *C ;
    Int *Super, *Map, *SuperMap ;
    size_t maxcsize ;
    Int nsuper, n, i, k, s, nrow ;
    int symbolic ;

    RETURN_IF_NULL_COMMON (FALSE) ;
    RETURN_IF_NULL (L, FALSE) ;
    RETURN_IF_NULL (A, FALSE) ;
    RETURN_IF_XTYPE_INVALID (A, CHOLMOD_REAL, CHOLMOD_ZOMPLEX, FALSE) ;
    RETURN_IF_XTYPE_INVALID (L, CHOLMOD_PATTERN, CHOLMOD_COMPLEX, FALSE) ;
    if (A->stype < 0)
    {
        if (A->nrow != A->ncol || A->nrow != L->n)
        {
            ERROR (CHOLMOD_INVALID, "invalid dimensions") ;
            return (FALSE) ;
        }
    }
    else if (A->stype == 0)
    {
        if (A->nrow != L->n)
        {
            ERROR (CHOLMOD_INVALID, "invalid dimensions") ;
            return (FALSE) ;
        }
        RETURN_IF_NULL (F, FALSE) ;
        RETURN_IF_XTYPE_INVALID (F, CHOLMOD_REAL, CHOLMOD_ZOMPLEX, FALSE) ;
        if (A->nrow != F->ncol || A->ncol != F->nrow || F->stype != 0)
        {
            ERROR (CHOLMOD_INVALID, "F invalid") ;
            return (FALSE) ;
        }
        if (A->xtype != F->xtype || A->dtype != F->dtype)
        {
            ERROR (CHOLMOD_INVALID, "A and F must have same xtype and dtype") ;
            return (FALSE) ;
        }
    }
    else
    {
        // symmetric upper case not supported
        ERROR (CHOLMOD_INVALID, "symmetric upper case not supported") ;
        return (FALSE) ;
    }
    if (!(L->is_super))
    {
        ERROR (CHOLMOD_INVALID, "L not supernodal") ;
        return (FALSE) ;
    }
    if (L->xtype != CHOLMOD_PATTERN)
    {
        if (! ((A->xtype == CHOLMOD_REAL    && L->xtype == CHOLMOD_REAL)
            || (A->xtype == CHOLMOD_COMPLEX && L->xtype == CHOLMOD_COMPLEX)
            || (A->xtype == CHOLMOD_ZOMPLEX && L->xtype == CHOLMOD_COMPLEX)))
        {
            ERROR (CHOLMOD_INVALID, "complex type mismatch") ;
            return (FALSE) ;
        }
        if (A->dtype != L->dtype)
        {
            ERROR (CHOLMOD_INVALID, "A and L must have the same dtype") ;
            return (FALSE) ;
        }
    }
    Common->status = CHOLMOD_OK ;

    //--------------------------------------------------------------------------
    // allocate workspace in Common
    //--------------------------------------------------------------------------

    nsuper = L->nsuper ;
    maxcsize = L->maxcsize ;
    nrow = A->nrow ;
    n = nrow ;

    PRINT1 (("nsuper "ID" maxcsize %g\n", nsuper, (double) maxcsize)) ;
    ASSERT (nsuper >= 0 && maxcsize > 0) ;

    // w = 2*nrow + 5*nsuper
    int ok = TRUE ;
    size_t w = CHOLMOD(mult_size_t) (A->nrow, 2, &ok) ;
    size_t t = CHOLMOD(mult_size_t) (L->nsuper, 5, &ok) ;
    w = CHOLMOD(add_size_t) (w, t, &ok) ;
    if (!ok)
    {
        ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
        return (FALSE) ;
    }

    CHOLMOD(allocate_work) (A->nrow, w, 0, Common) ;
    if (Common->status < CHOLMOD_OK)
    {
        return (FALSE) ;
    }
    ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 0, 0, Common)) ;

    //--------------------------------------------------------------------------
    // get the current factor L and allocate numerical part, if needed
    //--------------------------------------------------------------------------

    Super = L->super ;
    symbolic = (L->xtype == CHOLMOD_PATTERN) ;
    if (symbolic)
    {
        // convert to supernodal numeric by allocating L->x
        L->dtype = A->dtype ;       // ensure L has the same dtype as A
        CHOLMOD(change_factor) (
                (A->xtype == CHOLMOD_REAL) ? CHOLMOD_REAL : CHOLMOD_COMPLEX,
                TRUE, TRUE, TRUE, TRUE, L, Common) ;
        if (Common->status < CHOLMOD_OK)
        {
            // the factor L remains in symbolic supernodal form
            return (FALSE) ;
        }
    }
    ASSERT (L->dtype == A->dtype) ;
    ASSERT (L->dtype == CHOLMOD_DOUBLE || L->dtype == CHOLMOD_SINGLE) ;
    ASSERT (L->xtype == CHOLMOD_REAL   || L->xtype == CHOLMOD_COMPLEX) ;

    // supernodal LDL' is not supported
    L->is_ll = TRUE ;

    //--------------------------------------------------------------------------
    // get more workspace
    //--------------------------------------------------------------------------

    C = CHOLMOD(allocate_dense) (maxcsize, 1, maxcsize, L->xtype + L->dtype,
        Common) ;
    if (Common->status < CHOLMOD_OK)
    {
        int status = Common->status ;
        if (symbolic)
        {
            // Change L back to symbolic, since the numeric values are not
            // initialized.  This cannot fail.
            CHOLMOD(change_factor) (CHOLMOD_PATTERN, TRUE, TRUE, TRUE, TRUE, L,
                Common) ;
        }
        // the factor L is now back to the form it had on input
        Common->status = status ;
        return (FALSE) ;
    }

    //--------------------------------------------------------------------------
    // get workspace
    //--------------------------------------------------------------------------

    SuperMap = Common->Iwork ;  // size n
    Map = Common->Flag ;        // size n, use Flag as workspace for Map array
    for (i = 0 ; i < n ; i++)
    {
        Map [i] = EMPTY ;
    }

    //--------------------------------------------------------------------------
    // find the mapping of nodes to relaxed supernodes
    //--------------------------------------------------------------------------

    // SuperMap [k] = s if column k is contained in supernode s
    for (s = 0 ; s < nsuper ; s++)
    {
        PRINT1 (("Super ["ID"] "ID" ncols "ID"\n",
                    s, Super[s], Super[s+1]-Super[s])) ;
        for (k = Super [s] ; k < Super [s+1] ; k++)
        {
            SuperMap [k] = s ;
            PRINT2 (("relaxed SuperMap ["ID"] = "ID"\n", k, SuperMap [k])) ;
        }
    }

    //--------------------------------------------------------------------------
    // supernodal numerical factorization, using template routine
    //--------------------------------------------------------------------------

    float s_beta [2] ;
    s_beta [0] = (float) beta [0] ;
    s_beta [1] = (float) beta [1] ;

    switch ((A->xtype + A->dtype) % 8)
    {
        case CHOLMOD_REAL    + CHOLMOD_SINGLE:
            ok = rs_cholmod_super_numeric_worker (A, F, s_beta, L, C, Common) ;
            break ;

        case CHOLMOD_COMPLEX + CHOLMOD_SINGLE:
            ok = cs_cholmod_super_numeric_worker (A, F, s_beta, L, C, Common) ;
            break ;

        case CHOLMOD_ZOMPLEX + CHOLMOD_SINGLE:
            // A is zomplex, but L is complex
            ok = zs_cholmod_super_numeric_worker (A, F, s_beta, L, C, Common) ;
            break ;

        case CHOLMOD_REAL    + CHOLMOD_DOUBLE:
            ok = rd_cholmod_super_numeric_worker (A, F, beta, L, C, Common) ;
            break ;

        case CHOLMOD_COMPLEX + CHOLMOD_DOUBLE:
            ok = cd_cholmod_super_numeric_worker (A, F, beta, L, C, Common) ;
            break ;

        case CHOLMOD_ZOMPLEX + CHOLMOD_DOUBLE:
            // A is zomplex, but L is complex
            ok = zd_cholmod_super_numeric_worker (A, F, beta, L, C, Common) ;
            break ;
    }

    //--------------------------------------------------------------------------
    // clear Common workspace, free temp workspace C, and return
    //--------------------------------------------------------------------------

    // Flag array was used as workspace, clear it
    Common->mark = EMPTY ;
    CLEAR_FLAG (Common) ;
    ASSERT (check_flag (Common)) ;
    ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 0, 0, Common)) ;
    CHOLMOD(free_dense) (&C, Common) ;
    return (ok) ;
}

#endif
#endif

