/* ========================================================================== */
/* === Tcov/lpdemo ========================================================== */
/* ========================================================================== */
/* -----------------------------------------------------------------------------
* CHOLMOD/Tcov Module. Copyright (C) 2005-2006, Timothy A. Davis
* The CHOLMOD/Tcov 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
* -------------------------------------------------------------------------- */
/* A rectangular matrix is being tested (# nrows < # cols). This is a
* linear programming problem. Process the system using the same kind of
* operations that occur in an LP solver (the LP Dual Active Set Algorithm).
* This routine does not actually solve the LP. It simply mimics the kind
* of matrix operations that occur in LPDASA.
*
* The active set f is held in fset [0..fsize-1]. It is a subset of the columns
* of A. Columns not in the fset are in the list fnot [0..ncol-fsize-1].
*
* Rows can be added and deleted from A as well. A "dead" row is one that has
* been (temporarily) set to zero in A. If row i is dead, rflag [i] is 0,
* and 1 otherwise.
*
* The list r of "live" rows is kept in rset [0..rsize-1]. The list of "dead"
* rows is kept in rnot [0..nrow-rsize-1].
*
* The system to solve as r and/or f change is (beta*I + A(r,f)*A(r,f)') x = b.
* If a row i is deleted from A, it is set to zero. Row i of L and D are set
* to the ith row of the identity matrix.
*/
#include "cm.h"
#define MAXCOLS 8
/* ========================================================================== */
/* === Lcheck =============================================================== */
/* ========================================================================== */
/* Testing only: make sure there are no dead rows in L (excluding diagonal) */
static void Lcheck (cholmod_factor *L, Int *rflag)
{
Int *Lp, *Li, *Lnz ;
Int i, n, j, p, pend ;
double *Lx ;
if (L == NULL)
{
return ;
}
Lp = L->p ;
Li = L->i ;
Lx = L->x ;
Lnz = L->nz ;
n = L->n ;
for (j = 0 ; j < n ; j++)
{
p = Lp [j] ;
pend = p + Lnz [j] ;
for (p++ ; p < pend ; p++)
{
i = Li [p] ;
OK (IMPLIES (!rflag [i], Lx [p] == 0)) ;
}
}
}
/* ========================================================================== */
/* === lp_prune ============================================================= */
/* ========================================================================== */
/* C = A (r,f), except that C and A have the same row dimension. Row i of C
* and A(:,f) are equal if row i is in the rset. Row i of C is zero
* otherwise. C has as many columns as the size of f. */
cholmod_sparse *lp_prune
(
cholmod_sparse *A,
Int *rflag,
Int *fset,
Int fsize
)
{
cholmod_sparse *C ;
double *Ax, *Cx ;
Int *Ai, *Ap, *Ci, *Cp ;
Int i, kk, j, p, nz, nf, ncol ;
if (A == NULL)
{
ERROR (CHOLMOD_INVALID, "nothing to prune") ;
return (NULL) ;
}
Ap = A->p ;
Ai = A->i ;
Ax = A->x ;
ncol = A->ncol ;
nf = (fset == NULL) ? ncol : fsize ;
OK (fsize >= 0) ;
C = CHOLMOD(allocate_sparse) (A->nrow, nf, A->nzmax, A->sorted,
TRUE, 0, CHOLMOD_REAL, cm) ;
if (C == NULL)
{
ERROR (CHOLMOD_INVALID, "cannot create pruned C") ;
return (NULL) ;
}
Cp = C->p ;
Ci = C->i ;
Cx = C->x ;
nz = 0 ;
for (kk = 0 ; kk < nf ; kk++)
{
j = (fset == NULL) ? (kk) : (fset [kk]) ;
Cp [kk] = nz ;
for (p = Ap [j] ; p < Ap [j+1] ; p++)
{
i = Ai [p] ;
if (rflag [i])
{
Ci [nz] = i ;
Cx [nz] = Ax [p] ;
nz++ ;
}
}
}
Cp [nf] = nz ;
return (C) ;
}
/* ========================================================================== */
/* === lp_resid ============================================================= */
/* ========================================================================== */
/* Compute the 2-norm of the residual.
* norm ((beta*I + C*C')y(r) - b(r)), where C = A (r,f).
*/
double lp_resid
(
cholmod_sparse *A,
Int *rflag,
Int *fset,
Int fsize,
double beta [2],
cholmod_dense *Y,
cholmod_dense *B
)
{
cholmod_dense *R ;
double *Rx, *Yx ;
double rnorm, bnorm, ynorm, norm ;
cholmod_sparse *C ;
cholmod_dense *W ;
Int i, nrow ;
if (A == NULL)
{
ERROR (CHOLMOD_INVALID, "cannot compute LP resid") ;
return (1) ;
}
nrow = A->nrow ;
R = CHOLMOD(zeros) (nrow, 1, CHOLMOD_REAL, cm) ;
/* C = A(r,f). In LPDASA, we do this in place, without making a copy. */
C = lp_prune (A, rflag, fset, fsize) ;
/* W = C'*Y */
OK (fsize >= 0) ;
W = CHOLMOD(zeros) (fsize, 1, CHOLMOD_REAL, cm) ;
CHOLMOD(sdmult) (C, TRUE, one, zero, Y, W, cm) ;
/* R = B */
CHOLMOD(copy_dense2) (B, R, cm) ;
/* R = C*W - R */
CHOLMOD(sdmult) (C, FALSE, one, minusone, W, R, cm) ;
/* R = R + beta*Y, (beta = 1 for dropped rows) */
if (R != NULL && Y != NULL)
{
Rx = R->x ;
Yx = Y->x ;
for (i = 0 ; i < nrow ; i++)
{
if (rflag [i])
{
Rx [i] += beta [0] * Yx [i] ;
}
else
{
Rx [i] += Yx [i] ;
}
}
}
/* rnorm = norm (R) */
rnorm = CHOLMOD(norm_dense) (R, 2, cm) ;
bnorm = CHOLMOD(norm_dense) (B, 2, cm) ;
ynorm = CHOLMOD(norm_dense) (Y, 2, cm) ;
norm = MAX (bnorm, ynorm) ;
if (norm > 0)
{
rnorm /= norm ;
}
CHOLMOD(print_dense) (R, "R, resid", cm) ;
CHOLMOD(free_sparse) (&C, cm) ;
CHOLMOD(free_dense) (&W, cm) ;
CHOLMOD(free_dense) (&R, cm) ;
return (rnorm) ;
}
/* ========================================================================== */
/* === get_row ============================================================== */
/* ========================================================================== */
/* S = column i of beta*I + A(r,f)*A(r,f)' */
cholmod_sparse *get_row
(
cholmod_sparse *A,
Int i,
Int *rflag,
Int *fset,
Int fsize,
double beta [2]
)
{
cholmod_sparse *Ri, *R, *C, *S ;
double *Sx ;
Int *Sp, *Si ;
Int p, ii, found ;
if (rflag [i] == 0)
{
S = CHOLMOD(speye) (A->nrow, A->nrow, CHOLMOD_REAL, cm) ;
CHOLMOD(print_sparse) (S, "S identity", cm) ;
return (S) ;
}
OK (fsize >= 0) ;
/* Getting row i of A is expensive. In LPDASA, we maintain
* a copy of A(r,f)', and extact row i as column i of that
* matrix. We compute S = A(r,f)*A(i,f)' and S(i) += beta
* in a single pass. This is a simpler but slower method. */
/* R = A (i,f)' */
Ri = CHOLMOD(submatrix) (A, &i, 1, fset, fsize, TRUE, FALSE, cm) ;
R = CHOLMOD(transpose) (Ri, 1, cm) ;
CHOLMOD(free_sparse) (&Ri, cm) ;
/* C = A (r,f) */
C = lp_prune (A, rflag, fset, fsize) ;
/* S = C*R */
S = CHOLMOD(ssmult) (C, R, 0, TRUE, TRUE, cm) ;
CHOLMOD(free_sparse) (&C, cm) ;
CHOLMOD(free_sparse) (&R, cm) ;
if (S == NULL)
{
return (NULL) ;
}
/* S (i) += beta */
found = FALSE ;
Sp = S->p ;
Si = S->i ;
Sx = S->x ;
for (p = Sp [0] ; p < Sp [1] ; p++)
{
ii = Si [p] ;
if (ii == i)
{
found = TRUE ;
Sx [p] += beta [0] ;
break ;
}
}
if (!found)
{
/* oops, row index i is not present in S. Add it. */
CHOLMOD(reallocate_sparse) (S->nzmax+1, S, cm) ;
OK (Sp [1] < (Int) (S->nzmax)) ;
Si = S->i ;
Sx = S->x ;
Si [Sp [1]] = i ;
Sx [Sp [1]] = beta [0] ;
Sp [1]++ ;
S->sorted = FALSE ;
}
CHOLMOD(print_sparse) (S, "S", cm) ;
return (S) ;
}
/* ========================================================================== */
/* === lpdemo =============================================================== */
/* ========================================================================== */
double lpdemo (cholmod_triplet *T)
{
double r, maxerr = 0, anorm, bnorm, norm, xnorm, ynorm ;
double *b = NULL, *Yx = NULL, *Xx = NULL, *Sx ;
cholmod_sparse *A, *AT, *Apermuted, *C, *S, *Row ;
cholmod_dense *X, *B, *Y, *DeltaB, *R ;
cholmod_factor *L ;
Int *init, *rset, *rnot, *fset, *fnot, *rflag, *P, *Pinv, *Lperm, *fflag,
*Sp, *Si, *StaticParent ;
Int i, j, k, nrow, ncol, fsize, cols [MAXCOLS+1], trial, rank, kk, rsize,
p, op, ok ;
double beta [2], bk [2], yk [2] ;
/* ---------------------------------------------------------------------- */
/* convert T into a sparse matrix A */
/* ---------------------------------------------------------------------- */
if (T == NULL || T->ncol == 0)
{
/* nothing to do */
return (0) ;
}
if (T->xtype != CHOLMOD_REAL)
{
return (0) ;
}
A = CHOLMOD(triplet_to_sparse) (T, 0, cm) ;
if (A == NULL)
{
ERROR (CHOLMOD_INVALID, "cannot continue LP demo") ;
return (1) ;
}
nrow = A->nrow ;
ncol = A->ncol ;
anorm = CHOLMOD(norm_sparse) (A, 1, cm) ;
/* switch for afiro, but not galenet */
cm->supernodal_switch = 5 ;
/* ---------------------------------------------------------------------- */
/* select a random initial row and column basis */
/* ---------------------------------------------------------------------- */
/* select an initial fset of size nrow */
init = prand (ncol) ; /* RAND */
fset = CHOLMOD(malloc) (ncol, sizeof (Int), cm) ;
fnot = CHOLMOD(malloc) (ncol, sizeof (Int), cm) ;
fflag = CHOLMOD(malloc) (ncol, sizeof (Int), cm) ;
fsize = MIN (nrow, ncol) ;
if (init != NULL && fset != NULL && fflag != NULL)
{
for (k = 0 ; k < fsize ; k++)
{
j = init [k] ;
fset [k] = j ;
fflag [j] = 1 ;
}
for ( ; k < ncol ; k++)
{
j = init [k] ;
fnot [k-fsize] = j ;
fflag [j] = 0 ;
}
}
CHOLMOD(free) (ncol, sizeof (Int), init, cm) ;
/* all rows are live */
rsize = nrow ;
rflag = CHOLMOD(malloc) (nrow, sizeof (Int), cm) ;
rset = CHOLMOD(malloc) (nrow, sizeof (Int), cm) ;
rnot = CHOLMOD(malloc) (nrow, sizeof (Int), cm) ;
if (rset != NULL && rflag != NULL)
{
for (i = 0 ; i < nrow ; i++)
{
rflag [i] = 1 ;
rset [i] = i ;
}
}
/* ---------------------------------------------------------------------- */
/* factorize the first matrix, beta*I + A(p,f)*A(p,f)' */
/* ---------------------------------------------------------------------- */
beta [0] = 1e-6 ;
beta [1] = 0 ;
/* Need to prune entries due to relaxed amalgamation, or else
* cholmod_row_subtree will not be able to find all the entries in row
* k of L. */
cm->final_resymbol = TRUE ;
cm->final_asis = FALSE ;
cm->final_super = FALSE ;
cm->final_ll = FALSE ;
cm->final_pack = FALSE ;
cm->final_monotonic = FALSE ;
L = CHOLMOD(analyze_p) (A, NULL, fset, fsize, cm) ;
CHOLMOD(factorize_p) (A, beta, fset, fsize, L, cm) ;
/* get a copy of the fill-reducing permutation P and compute its inverse */
Lperm = (L != NULL) ? (L->Perm) : NULL ;
P = CHOLMOD(malloc) (nrow, sizeof (Int), cm) ;
Pinv = CHOLMOD(malloc) (nrow, sizeof (Int), cm) ;
if (P != NULL && Pinv != NULL && Lperm != NULL)
{
for (k = 0 ; k < nrow ; k++)
{
P [k] = Lperm [k] ;
Pinv [P [k]] = k ;
}
}
else
{
P = CHOLMOD(free) (nrow, sizeof (Int), P, cm) ;
Pinv = CHOLMOD(free) (nrow, sizeof (Int), Pinv, cm) ;
}
if (cm->print > 1)
{
k = cm->print ;
cm->print = 5 ;
CHOLMOD(print_common) ("cm for lpdemo", cm) ;
cm->print = k ;
}
/* ---------------------------------------------------------------------- */
/* A=P*A: permute the rows of A according to P */
/* ---------------------------------------------------------------------- */
/* This is done just once, since the system will be solved and modified
* many times. It's faster, and easier, to work in the permuted ordering
* rather than the original ordering. */
/* A will become unsorted later on; don't bother to sort it here */
Apermuted = CHOLMOD(submatrix) (A, P, nrow, NULL, -1, TRUE, TRUE, cm) ;
CHOLMOD(free_sparse) (&A, cm) ;
A = Apermuted ;
/* ---------------------------------------------------------------------- */
/* find the etree of A*A' */
/* ---------------------------------------------------------------------- */
/* Since the fset is a subset of 0:ncol-1, and rset is a subset of 0:nrow-1,
* the nonzero pattern of the Cholesky factorization of A(r,f)*A(r,f)' is a
* subset of the Cholesky factorization of A*A'. After many updates/
* downdates/rowadds/rowdels, any given row i of L may have entries that
* are not in the factorization of A (r,f)*A(r,f)'. To drop a row using
* cholmod_rowdel, we either need to know the pattern of the ith row of L,
* we can pass NULL and have cholmod_rowdel look at each column 0 to i-1.
* The StaticParent array is the etree of A*A', and it suffices to compute
* the pattern of the ith row of L based on that etree, and A and A'
* (ignoring the fset and rset). This gives us an upper bound on the
* nonzero pattern of the ith row of the current L (the factorization
* of A(r,f)*A(r,f)'.
*/
/* AT = nonzero pattern of A', used for row-subtree computations */
AT = CHOLMOD(transpose) (A, 0, cm) ;
/* Row = cholmod_row_subtree workspace (unsorted, packed, unsym, pattern) */
Row = CHOLMOD(allocate_sparse) (nrow, 1, nrow, FALSE, TRUE, 0,
CHOLMOD_PATTERN, cm) ;
/* Compute the "static" etree; the etree of A*A' */
StaticParent = CHOLMOD(malloc) (nrow, sizeof (Int), cm) ;
CHOLMOD(etree) (AT, StaticParent, cm) ;
/* ---------------------------------------------------------------------- */
/* compute initial right-hand-side */
/* ---------------------------------------------------------------------- */
/* If row i of the original A and B is row k of the permuted P*A and P*B,
* then P [k] = i and Pinv [i] = k. Row indices of A now refer to the
* permuted form of A, not the original A. Likewise, row k of B will
* refer to the permuted row k = Pinv [i], not the original row i. In a
* real program, this would affect how B is computed. This program just
* creates a random B anyway, so the order of B does not matter. It does
* use Pinv [i], just to show you how you would do it.
*/
B = CHOLMOD(zeros) (nrow, 1, CHOLMOD_REAL, cm) ;
if (B != NULL && Pinv != NULL)
{
b = B->x ;
for (i = 0 ; i < nrow ; i++)
{
/* row i of the original B is row k of the permuted B */
k = Pinv [i] ;
b [k] = xrand (1.) ; /* RAND */
}
}
/* ---------------------------------------------------------------------- */
/* solve the system */
/* ---------------------------------------------------------------------- */
/* Solve the system (beta*I + A(:,f)*A(:,f)')y=b without using L->Perm,
* since A and B have already been permuted according to L->Perm. */
DeltaB = CHOLMOD(zeros) (nrow, 1, CHOLMOD_REAL, cm) ;
/* solve Lx=b */
X = CHOLMOD(solve) (CHOLMOD_L, L, B, cm) ;
/* solve DL'y=x */
Y = CHOLMOD(solve) (CHOLMOD_DLt, L, X, cm) ;
r = lp_resid (A, rflag, fset, fsize, beta, Y, B) ;
MAXERR (maxerr, r, 1) ;
bk [0] = 0 ;
bk [1] = 0 ;
yk [0] = 0 ;
yk [1] = 0 ;
bnorm = CHOLMOD(norm_dense) (B, 1, cm) ;
/* ---------------------------------------------------------------------- */
/* modify the system */
/* ---------------------------------------------------------------------- */
ok = (fset != NULL && fnot != NULL && fflag != NULL &&
rset != NULL && rnot != NULL && rflag != NULL &&
B != NULL && Y != NULL && X != NULL && Row != NULL && A != NULL &&
AT != NULL && StaticParent != NULL && DeltaB != NULL && L != NULL &&
L->xtype != CHOLMOD_PATTERN && !(L->is_ll) && !(L->is_super)) ;
for (trial = 1 ; ok && trial < MAX (64, 2*ncol) ; trial++)
{
/* select an operation at random */
op = nrand (6) ; /* RAND */
Xx = X->x ;
Yx = Y->x ;
switch (op)
{
/* -------------------------------------------------------------- */
case 0: /* update */
/* -------------------------------------------------------------- */
/* pick some columns at random, but not all columns */
rank = 1 + nrand (MAXCOLS+4) ; /* RAND */
rank = MIN (rank, MAXCOLS) ;
rank = MIN (rank, ncol-fsize-1) ;
if (rank <= 0)
{
continue ;
}
/* remove the columns from fnot and add them to fset */
for (k = 0 ; k < rank ; k++)
{
kk = nrand (ncol-fsize) ; /* RAND */
j = fnot [kk] ;
fnot [kk] = fnot [ncol-fsize-1] ;
fset [fsize++] = j ;
OK (fsize < ncol) ;
cols [k] = j ;
fflag [j] = 1 ;
}
/* update L, and the solution to Lx=b+deltaB */
C = lp_prune (A, rflag, cols, rank) ;
ok = CHOLMOD(updown_solve) (TRUE, C, L, X, DeltaB, cm) ;
CHOLMOD(free_sparse) (&C, cm) ;
break ;
/* -------------------------------------------------------------- */
case 1: /* downdate */
/* -------------------------------------------------------------- */
/* pick some columns at random, but not all columns */
rank = 1 + nrand (MAXCOLS+4) ; /* RAND */
rank = MIN (rank, MAXCOLS) ;
rank = MIN (rank, fsize-1) ;
if (rank <= 0)
{
continue ;
}
/* remove the columns from fset and add them to fnot */
for (k = 0 ; k < rank ; k++)
{
kk = nrand (fsize) ; /* RAND */
j = fset [kk] ;
fset [kk] = fset [fsize-1] ;
fnot [ncol-fsize] = j ;
fsize-- ;
OK (fsize > 0) ;
cols [k] = j ;
fflag [j] = 0 ;
}
/* downdate L, and the solution to Lx=b+deltaB */
C = lp_prune (A, rflag, cols, rank) ;
ok = CHOLMOD(updown_solve) (FALSE, C, L, X, DeltaB, cm) ;
CHOLMOD(free_sparse) (&C, cm) ;
break ;
/* -------------------------------------------------------------- */
case 2: /* resymbol (no change to numerical values) */
/* -------------------------------------------------------------- */
/* let resymbol handle the fset */
C = lp_prune (A, rflag, NULL, 0) ;
ok = CHOLMOD(resymbol_noperm) (C, fset, fsize, TRUE, L, cm) ;
CHOLMOD(free_sparse) (&C, cm) ;
break;
/* -------------------------------------------------------------- */
case 3: /* add row */
/* -------------------------------------------------------------- */
/* remove a row from rnot and add to rset */
if (nrow == rsize)
{
continue ;
}
kk = nrand (nrow-rsize) ; /* RAND */
i = rnot [kk] ;
OK (rflag [i] == 0) ;
rnot [kk] = rnot [nrow-rsize-1] ;
rset [rsize++] = i ;
rflag [i] = 1 ;
/* S = column i of beta*I + A(r,f)*A(r,f)' */
S = get_row (A, i, rflag, fset, fsize, beta) ;
ok = (S != NULL) ;
if (ok)
{
/* pick a random right-hand-side for this new row */
b [i] = 1 ; /* xrand (1) */ /* was RAND */
bk [0] = b [i] ;
bk [1] = 0 ;
ok = CHOLMOD(rowadd_solve) (i, S, bk, L, X, DeltaB, cm) ;
}
CHOLMOD(free_sparse) (&S, cm) ;
break ;
/* -------------------------------------------------------------- */
case 4: /* delete row */
/* -------------------------------------------------------------- */
/* remove a row from rset and add to rnot */
if (rsize == 0)
{
continue ;
}
kk = nrand (rsize) ; /* RAND */
i = rset [kk] ;
OK (rflag [i] == 1) ;
rset [kk] = rset [rsize-1] ;
rnot [nrow-rsize] = i ;
rsize-- ;
/* S = column i of beta*I + A(r,f)*A(r,f)' */
S = get_row (A, i, rflag, fset, fsize, beta) ;
ok = (S != NULL) ;
if (ok)
{
/* B = B - S * y(i) */
Sp = S->p ;
Si = S->i ;
Sx = S->x ;
for (p = 0 ; p < Sp [1] ; p++)
{
b [Si [p]] -= Sx [p] * Yx [i] ;
}
/* B(i) = y(i) */
b [i] = Yx [i] ;
yk [0] = Yx [i] ;
yk [1] = 0 ;
/* pick a method arbitrarily */
if (trial % 2)
{
/* get upper bound nonzero pattern of L(i,0:i-1) */
CHOLMOD(row_subtree) (A, AT, i, StaticParent, Row, cm) ;
ok = CHOLMOD(rowdel_solve) (i, Row, yk, L, X, DeltaB,
cm) ;
}
else
{
/* Look in all cols 0 to i-1 for entries in L(i,0:i-1).
* This is more costly, but requires no knowledge of
* an upper bound on the pattern of L. */
ok = CHOLMOD(rowdel_solve) (i, NULL, yk, L, X, DeltaB,
cm) ;
}
/* for testing only, to ensure cholmod_row_subtree worked */
if (ok)
{
rflag [i] = 0 ;
Lcheck (L, rflag) ;
}
}
if (ok)
{
/* let resymbol handle the fset */
C = lp_prune (A, rflag, NULL, 0) ;
ok = CHOLMOD(resymbol_noperm) (C, fset, fsize, TRUE, L, cm) ;
CHOLMOD(free_sparse) (&C, cm) ;
}
CHOLMOD(free_sparse) (&S, cm) ;
break ;
/* -------------------------------------------------------------- */
case 5: /* convert, just for testing */
/* -------------------------------------------------------------- */
/* convert to LDL', optionally packed */
if (trial % 2)
{
ok = CHOLMOD(change_factor) (CHOLMOD_REAL, FALSE, FALSE,
TRUE, TRUE, L, cm) ;
}
else
{
ok = CHOLMOD(change_factor) (CHOLMOD_REAL, FALSE, FALSE,
FALSE, TRUE, L, cm) ;
}
break ;
}
if (ok)
{
/* scale B and X if their norm is getting large */
ynorm = CHOLMOD(norm_dense) (Y, 1, cm) ;
bnorm = CHOLMOD(norm_dense) (B, 1, cm) ;
xnorm = CHOLMOD(norm_dense) (X, 1, cm) ;
norm = MAX (bnorm, xnorm) ;
norm = MAX (norm, ynorm) ;
if (norm > 1e10)
{
for (i = 0 ; i < nrow ; i++)
{
Xx [i] /= norm ;
b [i] /= norm ;
}
}
CHOLMOD(free_dense) (&Y, cm) ;
Y = CHOLMOD(solve) (CHOLMOD_DLt, L, X, cm) ;
r = lp_resid (A, rflag, fset, fsize, beta, Y, B) ;
OK (!ISNAN (r)) ;
MAXERR (maxerr, r, 1) ;
if (r > 1e-6 && cm->print > 1)
{
printf ("lp err %.1g operation: "ID" ok "ID"\n", r, op, ok) ;
}
ok = (Y != NULL) ;
}
}
CHOLMOD(free_dense) (&Y, cm) ;
OK (CHOLMOD(print_common) ("cm in lpdemo", cm)) ;
/* ---------------------------------------------------------------------- */
/* convert to LDL packed, LDL unpacked or LL packed and solve again */
/* ---------------------------------------------------------------------- */
/* solve the new system and check the residual */
CHOLMOD(print_factor) (L, "L final, for convert", cm) ;
if (ok)
{
switch (nrand (3)) /* RAND */
{
/* pick one at random */
case 0:
{
ok = CHOLMOD(change_factor) (CHOLMOD_REAL, FALSE, FALSE, TRUE,
TRUE, L, cm) ;
Y = CHOLMOD(solve) (CHOLMOD_DLt, L, X, cm) ;
break ;
}
case 1:
{
ok = CHOLMOD(change_factor) (CHOLMOD_REAL, FALSE, FALSE, FALSE,
TRUE, L, cm) ;
Y = CHOLMOD(solve) (CHOLMOD_DLt, L, X, cm) ;
break ;
}
case 2:
{
ok = CHOLMOD(change_factor) (CHOLMOD_REAL, TRUE, FALSE, TRUE,
TRUE, L, cm) ;
Y = CHOLMOD(solve) (CHOLMOD_LDLt, L, B, cm) ;
break ;
}
}
r = lp_resid (A, rflag, fset, fsize, beta, Y, B) ;
OK (!ISNAN (r)) ;
MAXERR (maxerr, r, 1) ;
CHOLMOD(print_factor) (L, "L after convert", cm) ;
}
/* ---------------------------------------------------------------------- */
/* rank-1 update, but only partial Lx=b update */
/* ---------------------------------------------------------------------- */
if (ok && fsize < ncol && nrow > 3)
{
Int colmark [1] ;
j = fnot [0] ;
fnot [0] = fnot [ncol-fsize-1] ;
fset [fsize++] = j ;
OK (fsize <= ncol) ;
cols [0] = j ;
fflag [j] = 1 ;
for (colmark [0] = 0 ; colmark [0] <= nrow ; colmark [0]++)
{
cholmod_factor *L2 ;
cholmod_dense *X2 ;
double *X2x ;
L2 = CHOLMOD(copy_factor) (L, cm) ;
X2 = CHOLMOD(copy_dense) (X, cm) ;
X2x = (X2 == NULL) ? NULL : X2->x ;
/* fprintf (stderr, "check colmark "ID"\n", colmark [0]) ; */
printf ("check cholmark "ID"\n", colmark [0]) ;
/* colmark [0] = 3 ; */
/* update L, and the solution to Lx=b+deltaB,
* but only update solution in rows 0 to colmark[0] */
C = lp_prune (A, rflag, cols, 1) ;
ok = CHOLMOD(updown_mark) (TRUE, C, colmark, L2, X2, DeltaB, cm) ;
CHOLMOD(free_sparse) (&C, cm) ;
/* compare with Lr=b+deltaB */
R = CHOLMOD(solve) (CHOLMOD_L, L2, B, cm) ;
r = -1 ;
if (ok && R != NULL)
{
double *Rx ;
Rx = R->x ;
r = 0 ;
for (i = 0 ; i < colmark [0] ; i++)
{
r = MAX (r, fabs (X2x [i] - Rx [i])) ;
}
MAXERR (maxerr, r, 1) ;
}
printf ("check cholmark resid %6.2e\n", r) ;
CHOLMOD(free_dense) (&R, cm) ;
CHOLMOD(free_dense) (&X2, cm) ;
CHOLMOD(free_factor) (&L2, cm) ;
}
}
/* ---------------------------------------------------------------------- */
/* free everything */
/* ---------------------------------------------------------------------- */
/* restore defaults */
cm->final_resymbol = FALSE ;
cm->final_asis = TRUE ;
cm->supernodal_switch = 40 ;
CHOLMOD(free) (nrow, sizeof (Int), StaticParent, cm) ;
CHOLMOD(free) (nrow, sizeof (Int), Pinv, cm) ;
CHOLMOD(free) (nrow, sizeof (Int), P, cm) ;
CHOLMOD(free) (nrow, sizeof (Int), rflag, cm) ;
CHOLMOD(free) (nrow, sizeof (Int), rset, cm) ;
CHOLMOD(free) (nrow, sizeof (Int), rnot, cm) ;
CHOLMOD(free) (ncol, sizeof (Int), fset, cm) ;
CHOLMOD(free) (ncol, sizeof (Int), fnot, cm) ;
CHOLMOD(free) (ncol, sizeof (Int), fflag, cm) ;
CHOLMOD(free_factor) (&L, cm) ;
CHOLMOD(free_sparse) (&Row, cm) ;
CHOLMOD(free_sparse) (&AT, cm) ;
CHOLMOD(free_sparse) (&A, cm) ;
CHOLMOD(free_dense) (&B, cm) ;
CHOLMOD(free_dense) (&X, cm) ;
CHOLMOD(free_dense) (&Y, cm) ;
CHOLMOD(free_dense) (&DeltaB, cm) ;
progress (0, '.') ;
return (maxerr) ;
}
syntax highlighted by Code2HTML, v. 0.9.1