/* ========================================================================== */
/* === 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
syntax highlighted by Code2HTML, v. 0.9.1