//------------------------------------------------------------------------------
// CHOLMOD/Tcov/t_raw_factor: test CHOLMOD factorization and solvers
//------------------------------------------------------------------------------

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

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

// Factorize A using cholmod_rowfac for the simplicial case, and the
// cholmod_super_* routines for the supernodal case, and test the solution to
// linear systems.

#include "cm.h"

//------------------------------------------------------------------------------
// icomp
//------------------------------------------------------------------------------

// for sorting by qsort
static int icomp (Int *i, Int *j)
{
    if (*i < *j)
    {
        return (-1) ;
    }
    else
    {
        return (1) ;
    }
}

//------------------------------------------------------------------------------
// add_gunk
//------------------------------------------------------------------------------

static cholmod_sparse *add_gunk (cholmod_sparse *A)
{
    cholmod_sparse *S ;
    Real *Sx, *Sz ;
    Int *Sp, *Si, nz, p, save3, j, n ;

    if (A == NULL) return (NULL) ;

    A->nzmax++ ;
    S = CHOLMOD(copy_sparse) (A, cm) ;
    A->nzmax-- ;

    // add a S(n,1)=1 entry to the matrix
    if (S != NULL)
    {
        S->sorted = FALSE ;
        Sx = S->x ;
        Si = S->i ;
        Sp = S->p ;
        Sz = S->z ;
        n = S->ncol ;
        nz = Sp [n] ;
        for (j = 1 ; j <= n ; j++)
        {
            Sp [j]++ ;
        }
        if (S->xtype == CHOLMOD_REAL)
        {
            for (p = nz-1 ; p >= 0 ; p--)
            {
                Si [p+1] = Si [p] ;
                Sx [p+1] = Sx [p] ;
            }
            Si [0] = n-1 ;
            Sx [0] = 99999 ;
        }
        else if (S->xtype == CHOLMOD_COMPLEX)
        {
            for (p = nz-1 ; p >= 0 ; p--)
            {
                Si [p+1] = Si [p] ;
                Sx [2*p+2] = Sx [2*p] ;
                Sx [2*p+3] = Sx [2*p+1] ;
            }
            Si [0] = n-1 ;
            Sx [0] = 99999 ;
            Sx [1] = 0 ;
        }
        else if (S->xtype == CHOLMOD_ZOMPLEX)
        {
            for (p = nz-1 ; p >= 0 ; p--)
            {
                Si [p+1] = Si [p] ;
                Sx [p+1] = Sx [p] ;
                Sz [p+1] = Sz [p] ;
            }
            Si [0] = n-1 ;
            Sx [0] = 99999 ;
            Sz [0] = 0 ;
        }
    }

    return (S) ;
}

//------------------------------------------------------------------------------
// raw_factor
//------------------------------------------------------------------------------

// Factor A, without using any fill-reducing permutation.  This may fail due
// to catastrophic fill-in (which is the desired test result for a large
// arrowhead matrix).

double raw_factor (cholmod_sparse *A, Int check_errors)
{
    double maxerr = 0, r, anorm ;
    cholmod_sparse *AT, *C, *LT, *Lsparse, *S, *ST, *R, *A1 ;
    cholmod_factor *L, *Lcopy ;
    cholmod_dense *X, *W, *B, *X2 ;
    Int i, k, n, ok, ok1, ok2, trial, rnz, lnz, Lxtype, Axtype, posdef,
        prefer_zomplex, Bxtype ;
    Int *Parent, *Post, *First, *Level, *Ri, *Rp, *LTp = NULL, *LTi = NULL, *P,
        *mask, *RLinkUp ;
    int64_t lr ;
    double beta [2] ;
    uint64_t save ;

    //--------------------------------------------------------------------------
    // create the problem
    //--------------------------------------------------------------------------

    if (A == NULL || A->stype != 1)
    {
        return (0) ;
    }

    W = NULL ;
    X2 = NULL ;
    L = NULL ;
    n = A->nrow ;
    B = rhs (A, 1, n, 0) ;
    AT = CHOLMOD(transpose) (A, 2, cm) ;
    Parent = CHOLMOD(malloc) (n, sizeof (Int), cm) ;
    Post = CHOLMOD(malloc) (n, sizeof (Int), cm) ;
    First = CHOLMOD(malloc) (n, sizeof (Int), cm) ;
    Level = CHOLMOD(malloc) (n, sizeof (Int), cm) ;
    beta [0] = 0 ;
    beta [1] = 0 ;
    anorm = CHOLMOD(norm_sparse) (A, 1, cm) ;

    prefer_zomplex = (A->xtype == CHOLMOD_ZOMPLEX) ;
    Bxtype = A->xtype ;

    //--------------------------------------------------------------------------
    // supernodal factorization
    //--------------------------------------------------------------------------

    L = CHOLMOD(alloc_factor) (n, DTYPE, cm) ;
    ok1 = CHOLMOD(etree) (A, Parent, cm) ;
    lr = CHOLMOD(postorder) (Parent, n, NULL, Post, cm) ;
    ok2 = CHOLMOD(rowcolcounts) (AT, NULL, 0, Parent, Post,
        NULL, (L != NULL) ? (L->ColCount) : NULL, First, Level, cm) ;

    if (ok2)
    {
        printf ("raw_factor: cm->fl %g cm->lnz %g\n", cm->fl, cm->lnz) ;
    }

    if (check_errors)
    {
        OKP (AT) ;
        OKP (Parent) ;
        OKP (Post) ;
        OKP (First) ;
        OKP (Level) ;
        OK (AT->stype == -1) ;
        OKP (L) ;
        OK (ok1) ;
        OK (ok2) ;
        OK (lr >= 0) ;

        // rowcolcounts requires A in symmetric lower form
        ok = CHOLMOD(rowcolcounts) (A, NULL, 0, Parent, Post,
            NULL, L->ColCount, First, Level, cm) ;                  NOT (ok) ;
    }

    // super_symbolic needs A in upper form, so this will succeed
    // unless the problem is huge
    ok = CHOLMOD(super_symbolic) (A, NULL, Parent, L, cm) ;

    // super_symbolic should fail if lnz is too large
    if (cm->lnz > SIZE_MAX / 2)
    {
        printf ("raw_factor: problem is huge\n") ;
        NOT (ok) ;
        OK (L->xtype == CHOLMOD_PATTERN && !(L->is_super)) ;

        // try changing to LDL packed, which should also fail
        ok = CHOLMOD(change_factor) (CHOLMOD_REAL, FALSE, FALSE, TRUE, TRUE,
                L, cm) ;
        NOT (ok) ;

        CHOLMOD(free_factor) (&L, cm) ;
        CHOLMOD(free_sparse) (&AT, cm) ;
        CHOLMOD(free_dense) (&B, cm) ;
        CHOLMOD(free) (n, sizeof (Int), First, cm) ;
        CHOLMOD(free) (n, sizeof (Int), Level, cm) ;
        CHOLMOD(free) (n, sizeof (Int), Parent, cm) ;
        CHOLMOD(free) (n, sizeof (Int), Post, cm) ;
        return (0) ;
    }

    if (check_errors)
    {

        if (cm->status == CHOLMOD_OUT_OF_MEMORY
         || cm->status == CHOLMOD_TOO_LARGE)
        {
            // no test case will reach here, but check just to be safe
            printf ("raw_factor: out of memory for symbolic case %d\n",
                cm->status) ;
            CHOLMOD(free_factor) (&L, cm) ;
            CHOLMOD(free_sparse) (&AT, cm) ;
            CHOLMOD(free_dense) (&B, cm) ;
            CHOLMOD(free) (n, sizeof (Int), First, cm) ;
            CHOLMOD(free) (n, sizeof (Int), Level, cm) ;
            CHOLMOD(free) (n, sizeof (Int), Parent, cm) ;
            CHOLMOD(free) (n, sizeof (Int), Post, cm) ;
            return (0) ;
        }

        OK (ok) ;
        ok = CHOLMOD(super_symbolic)(A, NULL, Parent, L, cm) ;      NOT (ok) ;
        ok = CHOLMOD(super_symbolic)(AT, NULL, Parent, L, cm) ;     NOT (ok) ;
        ok = CHOLMOD(super_symbolic)(NULL, NULL, Parent, L, cm) ;   NOT (ok) ;
        ok = CHOLMOD(super_symbolic)(A, NULL, NULL, L, cm) ;        NOT (ok) ;
        ok = CHOLMOD(super_symbolic)(A, NULL, Parent, NULL, cm) ;   NOT (ok) ;
    }

    // super_numeric needs A in lower form, so this will succeed unless
    // the problem is huge
    ok = CHOLMOD(super_numeric) (AT, NULL, Zero, L, cm) ;

    if (check_errors)
    {

        if (cm->status == CHOLMOD_OUT_OF_MEMORY)
        {
            // For the 64-bit case, the Matrix/a1 problem will reach here
            printf ("raw_factor: out of memory for numeric case\n") ;
            CHOLMOD(free_factor) (&L, cm) ;
            CHOLMOD(free_sparse) (&AT, cm) ;
            CHOLMOD(free_dense) (&B, cm) ;
            CHOLMOD(free) (n, sizeof (Int), First, cm) ;
            CHOLMOD(free) (n, sizeof (Int), Level, cm) ;
            CHOLMOD(free) (n, sizeof (Int), Parent, cm) ;
            CHOLMOD(free) (n, sizeof (Int), Post, cm) ;
            return (0) ;
        }

        OK (ok) ;
        ok = CHOLMOD(super_numeric)(A, NULL, Zero, L, cm) ;         NOT (ok) ;
        ok = CHOLMOD(super_numeric)(NULL, NULL, Zero, L, cm) ;      NOT (ok) ;
        ok = CHOLMOD(super_numeric)(AT, NULL, Zero, NULL, cm) ;     NOT (ok) ;
    }

    // solve
    Lxtype = (L == NULL) ? CHOLMOD_REAL : L->xtype ;
    W = CHOLMOD(zeros) (n, 1, Lxtype + DTYPE, cm) ;
    X = CHOLMOD(copy_dense) (B, cm) ;
    if (Bxtype == CHOLMOD_ZOMPLEX)
    {
        CHOLMOD(dense_xtype) (CHOLMOD_COMPLEX + DTYPE, X, cm) ;
    }

    CHOLMOD(print_factor) (L, "L for super l/ltsolve", cm) ;
    CHOLMOD(print_dense) (W, "W", cm) ;
    CHOLMOD(print_dense) (X, "X", cm) ;

    ok1 = CHOLMOD(super_lsolve) (L, X, W, cm) ;
    CHOLMOD(print_dense) (X, "X", cm) ;

    ok2 = CHOLMOD(super_ltsolve) (L, X, W, cm) ;
    CHOLMOD(print_dense) (X, "X", cm) ;

    if (Bxtype == CHOLMOD_ZOMPLEX)
    {
        CHOLMOD(dense_xtype) (CHOLMOD_ZOMPLEX + DTYPE, X, cm) ;
    }

    r = resid (A, X, B) ;
    MAXERR (maxerr, r, 1) ;

    if (Bxtype == CHOLMOD_ZOMPLEX)
    {
        CHOLMOD(dense_xtype) (CHOLMOD_COMPLEX + DTYPE, X, cm) ;
    }

    if (check_errors)
    {
        OKP (W) ;
        OKP (X) ;
        OK (ok1) ;
        OK (ok2) ;
        ok = CHOLMOD(super_lsolve) (NULL, X, W, cm) ;               NOT (ok) ;
        ok = CHOLMOD(super_ltsolve) (NULL, X, W, cm) ;              NOT (ok) ;
        ok = CHOLMOD(super_lsolve) (L, NULL, W, cm) ;               NOT (ok) ;
        ok = CHOLMOD(super_ltsolve) (L, NULL, W, cm) ;              NOT (ok) ;
        ok = CHOLMOD(super_lsolve) (L, X, NULL, cm) ;               NOT (ok) ;
        ok = CHOLMOD(super_ltsolve) (L, X, NULL, cm) ;              NOT (ok) ;

        if (L != NULL && L->maxesize > 1)
        {
            // W is too small
            ok = CHOLMOD(free_dense) (&W, cm) ;                     OK (ok) ;
            W = CHOLMOD(zeros) (1, 1, Lxtype + DTYPE, cm) ;
            OKP (W) ;
            ok = CHOLMOD(super_lsolve) (L, X, W, cm) ;              NOT (ok) ;
            ok = CHOLMOD(super_ltsolve) (L, X, W, cm) ;             NOT (ok) ;
            ok = CHOLMOD(free_dense) (&W, cm) ;                     OK (ok) ;
            W = CHOLMOD(zeros) (n, 1, Lxtype + DTYPE, cm) ;
            OKP (W) ;
        }

        // X2 has the wrong dimensions
        X2 = CHOLMOD(zeros) (n+1, 1, Lxtype + DTYPE, cm) ;
        OKP (X2) ;
        ok = CHOLMOD(super_lsolve) (L, X2, W, cm) ;                 NOT (ok) ;
        ok = CHOLMOD(super_ltsolve) (L, X2, W, cm) ;                NOT (ok) ;
        CHOLMOD(free_dense) (&X2, cm) ;
    }

    CHOLMOD(free_dense) (&X, cm) ;

    // X2 is n-by-0, which is OK
    X2 = CHOLMOD(zeros) (n, 0, Lxtype + DTYPE, cm) ;
    ok1 = CHOLMOD(super_lsolve) (L, X2, W, cm) ;
    ok2 = CHOLMOD(super_ltsolve) (L, X2, W, cm) ;
    CHOLMOD(free_dense) (&W, cm) ;
    CHOLMOD(free_dense) (&X2, cm) ;

    if (check_errors)
    {
        OK (ok1) ;
        OK (ok2) ;
        test_memory_handler ( ) ;
        my_tries = 0 ;
        ok = CHOLMOD(super_symbolic) (A, NULL, Parent, L, cm) ;     NOT (ok) ;
        ok = CHOLMOD(super_numeric) (AT, NULL, Zero, L, cm) ;       NOT (ok) ;
        normal_memory_handler ( ) ;
        cm->error_handler = NULL ;
    }

    // R = space for result of row_subtree and row_lsubtree
    R = CHOLMOD(allocate_sparse)(n, 1, n, FALSE, TRUE, 0,
        CHOLMOD_PATTERN + DTYPE, cm) ;

    //--------------------------------------------------------------------------
    // erroneous factorization
    //--------------------------------------------------------------------------

    // cannot use rowfac or row_lsubtree on a supernodal factorization
    if (check_errors && n > 0)
    {
        ok = CHOLMOD(rowfac) (A, NULL, beta, 0, 0, L, cm) ;         NOT (ok) ;
        ok = CHOLMOD(row_lsubtree) (A, &i, 0, n-1, L, R, cm) ;      NOT (ok) ;
    }

    //--------------------------------------------------------------------------
    // convert to simplicial LDL'
    //--------------------------------------------------------------------------

    CHOLMOD(change_factor) (Lxtype, FALSE, FALSE, TRUE, TRUE, L, cm) ;

    // remove entries due to relaxed supernodal amalgamation
    CHOLMOD(resymbol) (A, NULL, 0, TRUE, L, cm) ;

    // refactorize a numeric factor
    posdef = 0 ;   // unknown
    if (A != NULL && A->stype >= 0)
    {
        if (A->stype > 0 && A->packed)
        {
            S = add_gunk (A) ;
            CHOLMOD(rowfac) (S, NULL, beta, 0, n, L, cm) ;
            if (S && S->xtype == CHOLMOD_COMPLEX)
            {
                CHOLMOD(sparse_xtype) (CHOLMOD_ZOMPLEX + DTYPE, S, cm) ;
            }
            ok = CHOLMOD(free_sparse) (&S, cm) ;
            OK (ok) ;
        }
        else
        {
            CHOLMOD(rowfac) (A, NULL, beta, 0, n, L, cm) ;
        }
        posdef = (cm->status == CHOLMOD_OK) ;
    }

    // convert to a sparse matrix, and transpose L
    Lcopy = CHOLMOD(copy_factor)(L, cm) ;
    Lsparse = CHOLMOD(factor_to_sparse) (L, cm) ;

    LT = CHOLMOD(transpose) (Lsparse, 0, cm) ;
    CHOLMOD(free_sparse) (&Lsparse, cm) ;

    if (LT != NULL)
    {
        LTp = LT->p ;
        LTi = LT->i ;
        OK (LT->packed) ;
    }

    // remove the unit diagonal of LT
    CHOLMOD(band_inplace) (1, n, -1, LT, cm) ;

    // ST = pattern of A(p,p)'
    P = (L == NULL) ? NULL : L->Perm ;
    ST = CHOLMOD(ptranspose) (A, 0, P, NULL, 0, cm) ;

    // S = pattern of A(p,p)
    S = CHOLMOD(transpose) (ST, 0, cm) ;
    ok = CHOLMOD(free_sparse) (&ST, cm) ;

    if (R != NULL && LT != NULL && posdef && A != NULL && A->stype >= 0
            && S != NULL && Lcopy != NULL)
    {
        LTp = LT->p ;
        LTi = LT->i ;

        Ri = R->i ;
        Rp = R->p ;

        save = my_seed ( ) ;                                    // RAND
        for (trial = 0 ; trial < 30 ; trial++)
        {
            // pick a row at random
            i = nrand (n) ;                                     // RAND

            // compute R = pattern of L(i,0:i-1), using row subtrees
            ok = CHOLMOD(row_subtree) (S, NULL, i, Parent, R, cm) ;
            if (!ok)
            {
                break ;
            }
            rnz = Rp [1] ;

            // sort R
            qsort (Ri, rnz, sizeof (Int),
                    (int (*) (const void *, const void *)) icomp) ;

            // compare with ith column of L transpose
            lnz = LTp [i+1] - LTp [i] ;
            ok = TRUE ;
            for (k = 0 ; k < MIN (rnz,lnz) ; k++)
            {
                ok = ok && (Ri [k] == LTi [LTp [i] + k]) ;
            }
            OK (ok) ;
            OK (rnz == lnz) ;

            // compute R = pattern of L(i,0:i-1), using row lsubtrees
            ok = CHOLMOD(row_lsubtree) (S, NULL, 0, i, Lcopy, R, cm) ;
            if (!ok)
            {
                break ;
            }
            rnz = Rp [1] ;

            // sort R
            qsort (Ri, rnz, sizeof (Int),
                    (int (*) (const void *, const void *)) icomp) ;

            // compare with ith column of L transpose
            lnz = LTp [i+1] - LTp [i] ;
            ok = TRUE ;
            for (k = 0 ; k < MIN (rnz,lnz) ; k++)
            {
                // printf ("%d vs %d\n", Ri [k], LTi [LTp [i] + k]) ;
                ok = ok && (Ri [k] == LTi [LTp [i] + k]) ;
            }
            OK (ok) ;
            OK (rnz == lnz) ;

            // L is symbolic, so cholmod_lsubtree will fail
            if (check_errors)
            {
                ok = CHOLMOD(row_lsubtree) (S, NULL, 0, i, L, R, cm) ;
                NOT (ok) ;
            }

        }
        my_srand (save) ;                                       // RAND
    }

    ok = CHOLMOD(free_factor) (&L, cm) ;                            OK (ok) ;
    ok = CHOLMOD(free_factor) (&Lcopy, cm) ;                        OK (ok) ;
    ok = CHOLMOD(free_sparse) (&LT, cm) ;                           OK (ok) ;
    ok = CHOLMOD(free_sparse) (&R, cm) ;                            OK (ok) ;
    ok = CHOLMOD(free_sparse) (&S, cm) ;                            OK (ok) ;

    //--------------------------------------------------------------------------
    // simplicial LDL' or LL' factorization with no analysis
    //--------------------------------------------------------------------------

    for (trial = 0 ; trial <= check_errors ; trial++)
    {

        // create a simplicial symbolic factor
        L = CHOLMOD(alloc_factor) (n, DTYPE, cm) ;
        ok = TRUE ;
        Axtype = (A == NULL) ? CHOLMOD_REAL : A->xtype ;

        if (check_errors)
        {
            OKP (L) ;
            if (trial == 0)
            {
                // convert to packed LDL' first, then unpacked
                ok = CHOLMOD(change_factor) (Axtype, FALSE, FALSE, TRUE,
                        TRUE, L, cm) ;
                OK (ok);
                ok = CHOLMOD(change_factor) (Axtype, FALSE, FALSE, FALSE,
                        TRUE, L, cm) ;
                OK (ok) ;
            }
            else if (trial == 1)
            {
                ok = CHOLMOD(rowfac)(NULL, NULL, beta, 0, 0, L, cm) ; NOT (ok) ;
                ok = CHOLMOD(rowfac)(A, NULL, beta, 0, 0, NULL, cm) ; NOT (ok) ;
                ok = CHOLMOD(rowfac)(AT, NULL, beta, 0, 0, L, cm) ;   NOT (ok) ;
                if (n > 1)
                {
                    A1 = CHOLMOD(allocate_sparse)(1, 1, 1, TRUE, TRUE, 1,
                        CHOLMOD_PATTERN + DTYPE, cm) ;
                    OKP (A1) ;
                    ok = CHOLMOD(rowfac)(A1, NULL, beta, 0, 0, L, cm); NOT (ok);
                    ok = CHOLMOD(free_sparse)(&A1, cm) ;                OK (ok);
                }
            }
            else
            {
                // convert to symbolic LL'
                ok = CHOLMOD(change_factor) (CHOLMOD_PATTERN, TRUE, FALSE, TRUE,
                        TRUE, L, cm) ;
                OK (ok) ;
                OK (L->is_ll) ;
            }
        }

        // factor
        CHOLMOD(print_factor) (L, "L for rowfac", cm) ;
        CHOLMOD(print_sparse) (A, "A for rowfac", cm) ;

        cm->dbound = 1e-15 ;
        cm->sbound = 1e-6 ;
        for (k = 0 ; ok && k < n ; k++)
        {
            if (!CHOLMOD(rowfac) (A, NULL, beta, k, k+1, L, cm))
            {
                ok = FALSE ;
            }
            if (cm->status == CHOLMOD_NOT_POSDEF)
            {
                // LL' factorization failed; subsequent rowfac's should fail
                k++ ;
                ok = CHOLMOD(rowfac) (A, NULL, beta, k, k+1, L, cm) ;
                NOT (ok) ;
                ok = TRUE ;
            }
        }
        cm->dbound = 0 ;
        cm->sbound = 0 ;

        if (check_errors)
        {
            OK (ok) ;
            ok = CHOLMOD(rowfac) (A, NULL, beta, n, n+1, L, cm) ;    NOT (ok) ;
            ok = TRUE ;
        }

        // solve
        if (ok)
        {
            cm->prefer_zomplex = prefer_zomplex ;
            X = CHOLMOD(solve) (CHOLMOD_A, L, B, cm) ;

            cm->prefer_zomplex = FALSE ;
            r = resid (A, X, B) ;
            MAXERR (maxerr, r, 1) ;
            CHOLMOD(free_dense) (&X, cm) ;
        }

        CHOLMOD(free_factor) (&L, cm) ;
    }

    //--------------------------------------------------------------------------
    // factor again with entries in the (ignored) lower part A
    //--------------------------------------------------------------------------

    if (A->packed)
    {
        L = CHOLMOD(alloc_factor) (n, DTYPE, cm) ;
        C = add_gunk (A) ;

        CHOLMOD(rowfac) (C, NULL, beta, 0, n, L, cm) ;

        X = CHOLMOD(solve) (CHOLMOD_A, L, B, cm) ;

        r = resid (A, X, B) ;
        MAXERR (maxerr, r, 1) ;
        CHOLMOD(free_sparse) (&C, cm) ;
        CHOLMOD(free_factor) (&L, cm) ;
        CHOLMOD(free_dense) (&X, cm) ;
    }

    //--------------------------------------------------------------------------
    // factor again using rowfac_mask (for LPDASA only)
    //--------------------------------------------------------------------------

    r = raw_factor2 (A, 0., 0) ;
    MAXERR (maxerr, r, 1) ;

    r = raw_factor2 (A, 1e-16, 0) ;
    MAXERR (maxerr, r, 1) ;

    //--------------------------------------------------------------------------
    // free the problem
    //--------------------------------------------------------------------------

    CHOLMOD(free_sparse) (&AT, cm) ;
    CHOLMOD(free_dense) (&B, cm) ;
    CHOLMOD(free) (n, sizeof (Int), First, cm) ;
    CHOLMOD(free) (n, sizeof (Int), Level, cm) ;
    CHOLMOD(free) (n, sizeof (Int), Parent, cm) ;
    CHOLMOD(free) (n, sizeof (Int), Post, cm) ;
    progress (0, '.') ;
    return (maxerr) ;
}

//------------------------------------------------------------------------------
// raw_factor2
//------------------------------------------------------------------------------

// A->stype can be 0 (lower), 1 (upper) or 0 (unsymmetric).  In the first two
// cases, Ax=b is solved.  In the third, A*A'x=b is solved.  No analysis and no
// fill-reducing ordering is used.  Both simplicial LL' and LDL' factorizations
// are used (testing rowfac_mask, for LPDASA only).

double raw_factor2 (cholmod_sparse *A, double alpha, int domask)
{
    Int n, i, prefer_zomplex, is_ll, xtype, sorted, axtype, stype ;
    Int *mask = NULL, *RLinkUp = NULL, nz = 0 ;
    Int *Cp = NULL, added_gunk ;
    double maxerr = 0, r = 0 ;
    cholmod_sparse *AT = NULL, *C = NULL, *CT = NULL, *CC = NULL, *C2 = NULL ;
    cholmod_factor *L = NULL ;
    cholmod_dense *B = NULL, *X = NULL ;
    double beta [2] ;

    if (A == NULL)
    {
        return (0) ;
    }
    n = A->nrow ;
    if (n > 1000)
    {
        printf ("\nSkipping rowfac, matrix too large\n") ;
        return (0) ;
    }
    axtype = A->xtype ;
    beta [0] = alpha ;
    beta [1] = 0 ;

    prefer_zomplex = (A->xtype == CHOLMOD_ZOMPLEX) ;
    AT = CHOLMOD(transpose) (A, 2, cm) ;

    // ensure C has stype of 0 or 1.  Do not prune any entries
    stype = A->stype ;
    if (stype >= 0)
    {
        A->stype = 0 ;
        C = CHOLMOD(copy_sparse) (A, cm) ;
        A->stype = stype ;
        if (C) C->stype = stype ;
        CT = AT ;

    }
    else
    {
        C = AT ;
        A->stype = 0 ;
        CT = CHOLMOD(copy_sparse) (A, cm) ;
        A->stype = stype ;
        if (CT) CT->stype = stype ;

        // only domask if C is symmetric and upper part stored
        domask = FALSE ;

    }

    mask = CHOLMOD(malloc) (n, sizeof (Int), cm) ;
    RLinkUp = CHOLMOD(malloc) (n, sizeof (Int), cm) ;

    if (C && cm->status == CHOLMOD_OK)
    {
        for (i = 0 ; i < n ; i++)
        {
            mask [i] = -1 ;
            RLinkUp [i] = i+1 ;
        }
    }
    else
    {
        domask = FALSE ;
    }

    if (C && !(C->packed) && !(C->sorted))
    {
        // do not do the unpacked or unsorted cases
        domask = FALSE ;
    }

    // make a copy of C and add some gunk if stype > 0
    added_gunk = (C && C->stype > 0) ;
    if (added_gunk)
    {
        C2 = add_gunk (C) ;
    }
    else
    {
        C2 = CHOLMOD(copy_sparse) (C, cm) ;
    }

    CC = CHOLMOD(copy_sparse) (C2, cm) ;

    if (CC && domask)
    {
        Int *Cp, *Ci, p ;
        Real *Cx, *Cz ;

        // this implicitly sets the first row/col of C to zero, except diag.
        mask [0] = 1 ;

        // CC = C2, and then set the first row/col to zero, except diagonal
        Cp = CC->p ;
        Ci = CC->i ;
        Cx = CC->x ;
        Cz = CC->z ;
        nz = Cp [n] ;
        switch (C->xtype)
        {
            case CHOLMOD_REAL:
                for (p = 1 ; p < nz ; p++)
                {
                    if (Ci [p] == 0) Cx [p] = 0 ;
                }
                break ;

            case CHOLMOD_COMPLEX:
                for (p = 1 ; p < nz ; p++)
                {
                    if (Ci [p] == 0)
                    {
                        Cx [2*p  ] = 0 ;
                        Cx [2*p+1] = 0 ;
                    }
                }
                break ;

            case CHOLMOD_ZOMPLEX:
                for (p = 1 ; p < nz ; p++)
                {
                    if (Ci [p] == 0)
                    {
                        Cx [p] = 0 ;
                        Cz [p] = 0 ;
                    }
                }
                break ;
        }
    }

    B = rhs (CC, 1, n, 0) ;

    for (sorted = 1 ; sorted >= 0 ; sorted--)
    {

        if (!sorted)
        {
            if (C2 && !added_gunk) C2->sorted = FALSE ;
            if (C) C->sorted = FALSE ;
            if (CT) CT->sorted = FALSE ;
        }

        for (is_ll = 0 ; is_ll <= 1 ; is_ll++)
        {
            for (xtype = 0 ; xtype <= 1 ; xtype++)
            {

                L = CHOLMOD(alloc_factor) (n, DTYPE, cm) ;
                if (L) L->is_ll = is_ll ;

                if (xtype)
                {
                    CHOLMOD (change_factor) (axtype, is_ll, 0, 0, 1, L, cm) ;
                }

                CHOLMOD(rowfac_mask) (sorted ? C : C2,
                    CT, beta, 0, n, mask, RLinkUp, L, cm) ;

                cm->prefer_zomplex = prefer_zomplex ;
                X = CHOLMOD(solve) (CHOLMOD_A, L, B, cm) ;
                cm->prefer_zomplex = FALSE ;

                r = resid (CC, X, B) ;
                MAXERR (maxerr, r, 1) ;

                printf ("rowfac mask: resid is %g\n", r) ;

                CHOLMOD(free_factor) (&L, cm) ;
                CHOLMOD(free_dense) (&X, cm) ;
            }
        }
    }

    CHOLMOD(free) (n, sizeof (Int), mask, cm) ;
    CHOLMOD(free) (n, sizeof (Int), RLinkUp, cm) ;

    CHOLMOD(free_sparse) (&C2, cm) ;
    CHOLMOD(free_sparse) (&CC, cm) ;
    CHOLMOD(free_sparse) (&CT, cm) ;
    CHOLMOD(free_sparse) (&C, cm) ;
    CHOLMOD(free_dense) (&B, cm) ;

    return (maxerr) ;
}