The R Project SVN R-packages

/* ========================================================================== */
/* === Modify/cholmod_rowdel ================================================ */
/* ========================================================================== */

/* -----------------------------------------------------------------------------
 * CHOLMOD/Modify Module.
 * Copyright (C) 2005-2006, Timothy A. Davis and William W. Hager.
 * The CHOLMOD/Modify Module is licensed under Version 2.0 of the GNU
 * General Public License.  See gpl.txt for a text of the license.
 * CHOLMOD is also available under other licenses; contact authors for details.
 * http://www.cise.ufl.edu/research/sparse
 * -------------------------------------------------------------------------- */

/* Deletes a row and column from an LDL' factorization.  The row and column k
 * is set to the kth row and column of the identity matrix.  Optionally
 * downdates the solution to Lx=b.
 *
 * workspace: Flag (nrow), Head (nrow+1), W (nrow*2), Iwork (2*nrow)
 *
 * Only real matrices are supported (exception: since only the pattern of R
 * is used, it can have any valid xtype).
 */

#ifndef NMODIFY

#include "cholmod_internal.h"
#include "cholmod_modify.h"


/* ========================================================================== */
/* === cholmod_rowdel ======================================================= */
/* ========================================================================== */

/* Sets the kth row and column of L to be the kth row and column of the identity
 * matrix, and updates L(k+1:n,k+1:n) accordingly.   To reduce the running time,
 * the caller can optionally provide the nonzero pattern (or an upper bound) of
 * kth row of L, as the sparse n-by-1 vector R.  Provide R as NULL if you want
 * CHOLMOD to determine this itself, which is easier for the caller, but takes
 * a little more time.
 */

int CHOLMOD(rowdel)
(
    /* ---- input ---- */
    size_t k,       /* row/column index to delete */
    cholmod_sparse *R,  /* NULL, or the nonzero pattern of kth row of L */
    /* ---- in/out --- */
    cholmod_factor *L,  /* factor to modify */
    /* --------------- */
    cholmod_common *Common
)
{
    double yk [2] ;
    yk [0] = 0. ;
    yk [1] = 0. ;
    return (CHOLMOD(rowdel_mark) (k, R, yk, NULL, L, NULL, NULL, Common)) ;
}


/* ========================================================================== */
/* === cholmod_rowdel_solve ================================================= */
/* ========================================================================== */

/* Does the same as cholmod_rowdel, but also downdates the solution to Lx=b.
 * When row/column k of A is "deleted" from the system A*y=b, this can induce
 * a change to x, in addition to changes arising when L and b are modified.
 * If this is the case, the kth entry of y is required as input (yk) */

int CHOLMOD(rowdel_solve)
(
    /* ---- input ---- */
    size_t k,       /* row/column index to delete */
    cholmod_sparse *R,  /* NULL, or the nonzero pattern of kth row of L */
    double yk [2],  /* kth entry in the solution to A*y=b */
    /* ---- in/out --- */
    cholmod_factor *L,  /* factor to modify */
    cholmod_dense *X,   /* solution to Lx=b (size n-by-1) */
    cholmod_dense *DeltaB,  /* change in b, zero on output */
    /* --------------- */
    cholmod_common *Common
)
{
    return (CHOLMOD(rowdel_mark) (k, R, yk, NULL, L, X, DeltaB, Common)) ;
}


/* ========================================================================== */
/* === cholmod_rowdel_mark ================================================== */
/* ========================================================================== */

/* Does the same as cholmod_rowdel_solve, except only part of L is used in
 * the update/downdate of the solution to Lx=b.  This routine is an "expert"
 * routine.  It is meant for use in LPDASA only.
 *
 * if R == NULL then columns 0:k-1 of L are searched for row k.  Otherwise, it
 * searches columns in the set defined by the pattern of the first column of R.
 * This is meant to be the pattern of row k of L (a superset of that pattern is
 * OK too).  R must be a permutation of a subset of 0:k-1.
 */

int CHOLMOD(rowdel_mark)
(
    /* ---- input ---- */
    size_t kdel,    /* row/column index to delete */
    cholmod_sparse *R,  /* NULL, or the nonzero pattern of kth row of L */
    double yk [2],  /* kth entry in the solution to A*y=b */
    Int *colmark,   /* Int array of size 1.  See cholmod_updown.c */
    /* ---- in/out --- */
    cholmod_factor *L,  /* factor to modify */
    cholmod_dense *X,   /* solution to Lx=b (size n-by-1) */
    cholmod_dense *DeltaB,  /* change in b, zero on output */
    /* --------------- */
    cholmod_common *Common
)
{
    double dk, sqrt_dk, xk, dj, fl ;
    double *Lx, *Cx, *W, *Xx, *Nx ;
    Int *Li, *Lp, *Lnz, *Ci, *Rj, *Rp, *Iwork ;
    cholmod_sparse *C, Cmatrix ;
    Int j, p, pend, kk, lnz, n, Cp [2], do_solve, do_update, left, k,
    right, middle, i, klast, given_row, rnz ;
    size_t s ;
    int ok = TRUE ;

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

    RETURN_IF_NULL_COMMON (FALSE) ;
    RETURN_IF_NULL (L, FALSE) ;
    RETURN_IF_XTYPE_INVALID (L, CHOLMOD_PATTERN, CHOLMOD_REAL, FALSE) ;
    n = L->n ;
    k = kdel ;
    if (kdel >= L->n || k < 0)
    {
    ERROR (CHOLMOD_INVALID, "k invalid") ;
    return (FALSE) ;
    }
    if (R == NULL)
    {
    Rj = NULL ;
    rnz = EMPTY ;
    }
    else
    {
    RETURN_IF_XTYPE_INVALID (R, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
    if (R->ncol != 1 || R->nrow != L->n)
    {
        ERROR (CHOLMOD_INVALID, "R invalid") ;
        return (FALSE) ;
    }
    Rj = R->i ;
    Rp = R->p ;
    rnz = Rp [1] ;
    }
    do_solve = (X != NULL) && (DeltaB != NULL) ;
    if (do_solve)
    {
    RETURN_IF_XTYPE_INVALID (X, CHOLMOD_REAL, CHOLMOD_REAL, FALSE) ;
    RETURN_IF_XTYPE_INVALID (DeltaB, CHOLMOD_REAL, CHOLMOD_REAL, FALSE) ;
    Xx = X->x ;
    Nx = DeltaB->x ;
    if (X->nrow != L->n || X->ncol != 1 || DeltaB->nrow != L->n ||
        DeltaB->ncol != 1 || Xx == NULL || Nx == NULL)
    {
        ERROR (CHOLMOD_INVALID, "X and/or DeltaB invalid") ;
        return (FALSE) ;
    }
    }
    else
    {
    Xx = NULL ;
    Nx = NULL ;
    }
    Common->status = CHOLMOD_OK ;

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

    /* s = 2*n */
    s = CHOLMOD(mult_size_t) (n, 2, &ok) ;
    if (!ok)
    {
    ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
    return (FALSE) ;
    }

    CHOLMOD(allocate_work) (n, s, s, Common) ;
    if (Common->status < CHOLMOD_OK)
    {
    return (FALSE) ;
    }
    ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 2*n, Common)) ;

    /* ---------------------------------------------------------------------- */
    /* convert to simplicial numeric LDL' factor, if not already */
    /* ---------------------------------------------------------------------- */

    if (L->xtype == CHOLMOD_PATTERN || L->is_super || L->is_ll) 
    {
    /* can only update/downdate a simplicial LDL' factorization */
    CHOLMOD(change_factor) (CHOLMOD_REAL, FALSE, FALSE, FALSE, FALSE, L,
        Common) ;
    if (Common->status < CHOLMOD_OK)
    {
        /* out of memory, L is returned unchanged */
        return (FALSE) ;
    }
    }

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

    /* inputs, not modified on output: */
    Lp = L->p ;     /* size n+1 */

    /* outputs, contents defined on input for incremental case only: */
    Lnz = L->nz ;   /* size n */
    Li = L->i ;     /* size L->nzmax.  Can change in size. */
    Lx = L->x ;     /* size L->nzmax.  Can change in size. */

    ASSERT (L->nz != NULL) ;

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

    W = Common->Xwork ;     /* size n, used only in cholmod_updown */
    Cx = W + n ;        /* use 2nd column of Xwork for C (size n) */
    Iwork = Common->Iwork ;
    Ci = Iwork + n ;        /* size n (i/i/l) */
    /* NOTE: cholmod_updown uses Iwork [0..n-1] (i/i/l) as Stack */

    /* ---------------------------------------------------------------------- */
    /* prune row k from all columns of L */
    /* ---------------------------------------------------------------------- */

    given_row = (rnz >= 0) ;
    klast = given_row ? rnz : k ;
    PRINT2 (("given_row "ID"\n", given_row)) ;

    for (kk = 0 ; kk < klast ; kk++)
    {
    /* either search j = 0:k-1 or j = Rj [0:rnz-1] */
    j = given_row ? (Rj [kk]) : (kk) ;

    if (j < 0 || j >= k)
    {
        ERROR (CHOLMOD_INVALID, "R invalid") ;
        return (FALSE) ;
    }

    PRINT2 (("Prune col j = "ID":\n", j)) ;

    lnz = Lnz [j] ;
    dj = Lx [Lp [j]] ;
    ASSERT (Lnz [j] > 0 && Li [Lp [j]] == j) ;

    if (lnz > 1)
    {
        left = Lp [j] ;
        pend = left + lnz ;
        right = pend - 1 ;

        i = Li [right] ;

        if (i < k)
        {
        /* row k is not in column j */
        continue ;
        }
        else if (i == k)
        {
        /* k is the last row index in this column (quick delete) */
        if (do_solve)
        {
            Xx [j] -= yk [0] * dj * Lx [right] ;
        }
        Lx [right] = 0 ;
        }
        else
        {
        /* binary search for row k in column j */
        PRINT2 (("\nBinary search: lnz "ID" k = "ID"\n", lnz, k)) ;
        while (left < right)
        {
            middle = (left + right) / 2 ;
            PRINT2 (("left "ID" right "ID" middle "ID": ["ID" "ID""
            ""ID"]\n", left, right, middle,
            Li [left], Li [middle], Li [right])) ;
            if (k > Li [middle])
            {
            left = middle + 1 ;
            }
            else
            {
            right = middle ;
            }
        }
        ASSERT (left >= Lp [j] && left < pend) ;

#ifndef NDEBUG
        /* brute force, linear-time search */
        {
            Int p3 = Lp [j] ;
            i = EMPTY ;
            PRINT2 (("Brute force:\n")) ;
            for ( ; p3 < pend ; p3++)
            {
            i = Li [p3] ;
            PRINT2 (("p "ID" ["ID"]\n", p3, i)) ;
            if (i >= k)
            {
                break ;
            }
            }
            if (i == k)
            {
            ASSERT (k == Li [p3]) ;
            ASSERT (p3 == left) ;
            }
        }
#endif

        if (k == Li [left])
        {
            if (do_solve)
            {
            Xx [j] -= yk [0] * dj * Lx [left] ;
            }
            /* found row k in column j.  Prune it from the column.*/
            Lx [left] = 0 ;
        }
        }
    }
    }

#ifndef NDEBUG
    /* ensure that row k has been deleted from the matrix L */
    for (j = 0 ; j < k ; j++)
    {
    Int lasti ;
    lasti = EMPTY ;
    p = Lp [j] ;
    pend = p + Lnz [j] ;
    /* look for row k in column j */
    PRINT1 (("Pruned column "ID"\n", j)) ;
    for ( ; p < pend ; p++)
    {
        i = Li [p] ;
        PRINT2 ((" "ID"", i)) ;
        PRINT2 ((" %g\n", Lx [p])) ;
        ASSERT (IMPLIES (i == k, Lx [p] == 0)) ;
        ASSERT (i > lasti) ;
        lasti = i ;
    }
    PRINT1 (("\n")) ;
    }
#endif

    /* ---------------------------------------------------------------------- */
    /* set diagonal and clear column k of L */
    /* ---------------------------------------------------------------------- */

    lnz = Lnz [k] - 1 ;
    ASSERT (Lnz [k] > 0) ;

    /* ---------------------------------------------------------------------- */
    /* update/downdate */
    /* ---------------------------------------------------------------------- */

    /* update or downdate L (k+1:n, k+1:n) with the vector
     * C = L (:,k) * sqrt (abs (D [k]))
     * Do a numeric update if D[k] > 0, numeric downdate otherwise.
     */

    PRINT1 (("rowdel downdate lnz = "ID"\n", lnz)) ;

    /* store the new unit diagonal */
    p = Lp [k] ;
    pend = p + lnz + 1 ;
    dk = Lx [p] ;
    Lx [p++] = 1 ;
    PRINT2 (("D [k = "ID"] = %g\n", k, dk)) ;
    ok = TRUE ;
    fl = 0 ;

    if (lnz > 0)
    {
    /* compute DeltaB for updown (in DeltaB) */
    if (do_solve)
    {
        xk = Xx [k] - yk [0] * dk ;
        for ( ; p < pend ; p++)
        {
        Nx [Li [p]] += Lx [p] * xk ;
        }
    }

    do_update = IS_GT_ZERO (dk) ;
    if (!do_update)
    {
        dk = -dk ;
    }
    sqrt_dk = sqrt (dk) ;
    p = Lp [k] + 1 ;
    for (kk = 0 ; kk < lnz ; kk++, p++)
    {
        Ci [kk] = Li [p] ;
        Cx [kk] = Lx [p] * sqrt_dk ;
        Lx [p] = 0 ;        /* clear column k */
    }
    fl = lnz + 1 ;

    /* create a n-by-1 sparse matrix to hold the single column */
    C = &Cmatrix ;
    C->nrow = n ;
    C->ncol = 1 ;
    C->nzmax = lnz ;
    C->sorted = TRUE ;
    C->packed = TRUE ;
    C->p = Cp ;
    C->i = Ci ;
    C->x = Cx ;
    C->nz = NULL ;
    C->itype = L->itype ;
    C->xtype = L->xtype ;
    C->dtype = L->dtype ;
    C->z = NULL ;
    C->stype = 0 ;

    Cp [0] = 0 ;
    Cp [1] = lnz ;

    /* numeric update if dk > 0, and with Lx=b change */
    /* workspace: Flag (nrow), Head (nrow+1), W (nrow), Iwork (2*nrow) */
    ok = CHOLMOD(updown_mark) (do_update ? (1) : (0), C, colmark,
        L, X, DeltaB, Common) ;

    /* clear workspace */
    for (kk = 0 ; kk < lnz ; kk++)
    {
        Cx [kk] = 0 ;
    }
    }

    Common->modfl += fl ;

    if (do_solve)
    {
    /* kth equation becomes identity, so X(k) is now Y(k) */
    Xx [k] = yk [0] ;
    }

    DEBUG (CHOLMOD(dump_factor) (L, "LDL factorization, L:", Common)) ;
    ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 2*n, Common)) ;
    return (ok) ;
}
#endif