1 #include "cs_mex.h"
2 /* cs_updown: sparse Cholesky update/downdate (rank-1 or multiple rank) */
mexFunction(int nargout,mxArray * pargout[],int nargin,const mxArray * pargin[])3 void mexFunction
4 (
5 int nargout,
6 mxArray *pargout [ ],
7 int nargin,
8 const mxArray *pargin [ ]
9 )
10 {
11 cs Lmatrix, *Lin, Cmatrix, *C, *L, Cvector, *Cvec ;
12 csi ignore, j, k, n, lnz, *parent, sigma = 1, cp [2] ;
13 char sigma_string [20] ;
14
15 if (nargout > 1 || nargin < 3 || nargin > 4)
16 {
17 mexErrMsgTxt ("Usage: L = cs_updown(L,C,parent,sigma)") ;
18 }
19 Lin = cs_mex_get_sparse (&Lmatrix, 1, 1, pargin [0]) ; /* get input L */
20 n = Lin->n ;
21 if (nargin > 3 && mxIsChar (pargin [3]))
22 {
23 mxGetString (pargin [3], sigma_string, 8) ;
24 sigma = (sigma_string [0] == '-') ? (-1) : 1 ;
25 }
26 /* make a copy of L (this can take more work than updating L itself) */
27 lnz = Lin->p [n] ;
28 L = cs_spalloc (n, n, lnz, 1, 0) ;
29 for (j = 0 ; j <= n ; j++) L->p [j] = Lin->p [j] ;
30 for (k = 0 ; k < lnz ; k++) L->i [k] = Lin->i [k] ;
31 for (k = 0 ; k < lnz ; k++) L->x [k] = Lin->x [k] ;
32 cs_mex_check (0, n, -1, 0, 1, 1, pargin [1]) ; /* get C */
33 C = cs_mex_get_sparse (&Cmatrix, 0, 1, pargin [1]) ;
34 parent = cs_mex_get_int (n, pargin [2], &ignore, 0) ; /* get parent */
35
36 /* do the update one column at a time */
37 Cvec = &Cvector ;
38 Cvec->m = n ;
39 Cvec->n = 1 ;
40 Cvec->p = cp ;
41 Cvec->nz = -1 ;
42 cp [0] = 0 ;
43 for (k = 0 ; k < C->n ; k++)
44 {
45 /* extract C(:,k) */
46 cp [1] = C->p [k+1] - C->p [k] ;
47 Cvec->nzmax = cp [1] ;
48 Cvec->i = C->i + C->p [k] ;
49 Cvec->x = C->x + C->p [k] ;
50 cs_updown (L, sigma, Cvec, parent) ; /* update/downdate */
51 }
52 pargout [0] = cs_mex_put_sparse (&L) ; /* return new L */
53 }
54