//------------------------------------------------------------------------------
// CHOLMOD/Cholesky/cholmod_etree: elimination tree of A or A'*A
//------------------------------------------------------------------------------

// CHOLMOD/Cholesky Module.  Copyright (C) 2005-2023, Timothy A. Davis
// All Rights Reserved.
// SPDX-License-Identifier: LGPL-2.1+

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

// Compute the elimination tree of A or A'*A
//
// In the symmetric case, the upper triangular part of A is used.  Entries not
// in this part of the matrix are ignored.  Computing the etree of a symmetric
// matrix from just its lower triangular entries is not supported.
//
// In the unsymmetric case, all of A is used, and the etree of A'*A is computed.
//
// References:
//
// J. Liu, "A compact row storage scheme for Cholesky factors", ACM Trans.
// Math. Software, vol 12, 1986, pp. 127-148.
//
// J. Liu, "The role of elimination trees in sparse factorization", SIAM J.
// Matrix Analysis & Applic., vol 11, 1990, pp. 134-172.
//
// J. Gilbert, X. Li, E. Ng, B. Peyton, "Computing row and column counts for
// sparse QR and LU factorization", BIT, vol 41, 2001, pp. 693-710.
//
// workspace: symmetric: Iwork (nrow), unsymmetric: Iwork (nrow+ncol)
//
// Supports any xtype (pattern, real, complex, or zomplex) and any dtype.

#include "cholmod_internal.h"

#ifndef NCHOLESKY

//------------------------------------------------------------------------------
// update_etree
//------------------------------------------------------------------------------

static void update_etree
(
    // input:
    Int k,              // process the edge (k,i) in the input graph
    Int i,
    // input/output:
    Int Parent [ ],     // Parent [t] = p if p is the parent of t
    Int Ancestor [ ]    // Ancestor [t] is the ancestor of node t in the
                        // partially-constructed etree
)
{
    Int a ;
    for ( ; ; )     // traverse the path from k to the root of the tree
    {
        a = Ancestor [k] ;
        if (a == i)
        {
            // final ancestor reached; no change to tree
            return ;
        }
        // perform path compression
        Ancestor [k] = i ;
        if (a == EMPTY)
        {
            // final ancestor undefined; this is a new edge in the tree
            Parent [k] = i ;
            return ;
        }
        // traverse up to the ancestor of k
        k = a ;
    }
}

//------------------------------------------------------------------------------
// cholmod_etree
//------------------------------------------------------------------------------

// Find the elimination tree of A or A'*A

int CHOLMOD(etree)
(
    // input:
    cholmod_sparse *A,
    // output:
    Int *Parent,        // size ncol.  Parent [j] = p if p is the parent of j
    cholmod_common *Common
)
{

    Int *Ap, *Ai, *Anz, *Ancestor, *Prev, *Iwork ;
    Int i, j, jprev, p, pend, nrow, ncol, packed, stype ;

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

    RETURN_IF_NULL_COMMON (FALSE) ;
    RETURN_IF_NULL (A, FALSE) ;
    RETURN_IF_NULL (Parent, FALSE) ;
    RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
    Common->status = CHOLMOD_OK ;

    //--------------------------------------------------------------------------
    // allocate workspace
    //--------------------------------------------------------------------------

    stype = A->stype ;

    // s = A->nrow + (stype ? 0 : A->ncol)
    int ok = TRUE ;
    size_t s = CHOLMOD(add_size_t) (A->nrow, (stype ? 0 : A->ncol), &ok) ;
    if (!ok)
    {
        ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
        return (FALSE) ;
    }

    CHOLMOD(allocate_work) (0, s, 0, Common) ;
    if (Common->status < CHOLMOD_OK)
    {
        return (FALSE) ;        // out of memory
    }

    ASSERT (CHOLMOD(dump_sparse) (A, "etree", Common) >= 0) ;
    Iwork = Common->Iwork ;

    //--------------------------------------------------------------------------
    // get inputs
    //--------------------------------------------------------------------------

    ncol = A->ncol ;    // the number of columns of A
    nrow = A->nrow ;    // the number of rows of A
    Ap = A->p ;         // size ncol+1, column pointers for A
    Ai = A->i ;         // the row indices of A
    Anz = A->nz ;       // number of nonzeros in each column of A
    packed = A->packed ;
    Ancestor = Iwork ;  // size ncol

    for (j = 0 ; j < ncol ; j++)
    {
        Parent [j] = EMPTY ;
        Ancestor [j] = EMPTY ;
    }

    //--------------------------------------------------------------------------
    // compute the etree
    //--------------------------------------------------------------------------

    if (stype > 0)
    {

        //----------------------------------------------------------------------
        // symmetric (upper) case: compute etree (A)
        //----------------------------------------------------------------------

        for (j = 0 ; j < ncol ; j++)
        {
            // for each row i in column j of triu(A), excluding the diagonal
            p = Ap [j] ;
            pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
            for ( ; p < pend ; p++)
            {
                i = Ai [p] ;
                if (i < j)
                {
                    update_etree (i, j, Parent, Ancestor) ;
                }
            }
        }

    }
    else if (stype == 0)
    {

        //----------------------------------------------------------------------
        // unsymmetric case: compute etree (A'*A)
        //----------------------------------------------------------------------

        Prev = Iwork + ncol ;   // size nrow
        for (i = 0 ; i < nrow ; i++)
        {
            Prev [i] = EMPTY ;
        }
        for (j = 0 ; j < ncol ; j++)
        {
            // for each row i in column j of A
            p = Ap [j] ;
            pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
            for ( ; p < pend ; p++)
            {
                // a graph is constructed dynamically with one path per row
                // of A.  If the ith row of A contains column indices
                // (j1,j2,j3,j4) then the new graph has edges (j1,j2), (j2,j3),
                // and (j3,j4).  When at node i of this path-graph, all edges
                // (jprev,j) are considered, where jprev<j
                i = Ai [p] ;
                jprev = Prev [i] ;
                if (jprev != EMPTY)
                {
                    update_etree (jprev, j, Parent, Ancestor) ;
                }
                Prev [i] = j ;
            }
        }

    }
    else
    {

        //----------------------------------------------------------------------
        // symmetric case with lower triangular part not supported
        //----------------------------------------------------------------------

        ERROR (CHOLMOD_INVALID, "symmetric lower not supported") ;
        return (FALSE) ;
    }

    ASSERT (CHOLMOD(dump_parent) (Parent, ncol, "Parent", Common)) ;
    return (TRUE) ;
}
#endif