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_INT ignore, j, k, n, lnz, *parent, sigma = 1, cp [2], ok ;
12     char sigma_string [20] ;
13     if (nargout > 1 || nargin < 3 || nargin > 4)
14     {
15         mexErrMsgTxt ("Usage: L = cs_updown(L,C,parent,sigma)") ;
16     }
17     if (nargin > 3 && mxIsChar (pargin [3]))
18     {
19         mxGetString (pargin [3], sigma_string, 8) ;
20         sigma = (sigma_string [0] == '-') ? (-1) : 1 ;
21     }
22     n = mxGetN (pargin [0]) ;
23     parent = cs_dl_mex_get_int (n, pargin [2], &ignore, 0) ; /* get parent*/
24 
25     if (mxIsComplex (pargin [0]) || mxIsComplex (pargin [1]))
26     {
27 
28 #ifndef NCOMPLEX
29         cs_cl Lmatrix, *Lin, Cmatrix, *C, *L, Cvector, *Cvec ;
30         /* get input L, and copy MATLAB complex to C complex type */
31         Lin = cs_cl_mex_get_sparse (&Lmatrix, 1, pargin [0]) ;
32 
33         /* make a copy of L (this can take more work than updating L itself) */
34         lnz = Lin->p [n] ;
35         L = cs_cl_spalloc (n, n, lnz, 0, 0) ;
36         for (j = 0 ; j <= n ; j++) L->p [j] = Lin->p [j] ;
37         for (k = 0 ; k < lnz ; k++) L->i [k] = Lin->i [k] ;
38 
39         /* complex values already copied into Lin->x, use shallow copy for L */
40         L->x = Lin->x ;
41 
42         cs_mex_check (0, n, -1, 0, 1, 1, pargin [1]) ;     /* get C */
43         C = cs_cl_mex_get_sparse (&Cmatrix, 0, pargin [1]) ;
44 
45         /* do the update one column at a time */
46         Cvec = &Cvector ;
47         Cvec->m = n ;
48         Cvec->n = 1 ;
49         Cvec->p = cp ;
50         Cvec->nz = -1 ;
51         cp [0] = 0 ;
52         for (k = 0 ; k < C->n ; k++)
53         {
54             /* extract C(:,k) */
55             cp [1] = C->p [k+1] - C->p [k] ;
56             Cvec->nzmax = cp [1] ;
57             Cvec->i = C->i + C->p [k] ;
58             Cvec->x = C->x + C->p [k] ;
59             /* update/downdate */
60             ok = cs_cl_updown (L, sigma, Cvec, parent) ;
61             if (!ok) mexErrMsgTxt ("matrix is not positive definite") ;
62         }
63         /* return new L */
64         pargout [0] = cs_cl_mex_put_sparse (&L) ;
65 
66         cs_free (C->x) ;        /* free complex copy of C */
67 #else
68         mexErrMsgTxt ("complex matrices not supported") ;
69 #endif
70 
71     }
72     else
73     {
74 
75         cs_dl Lmatrix, *Lin, Cmatrix, *C, *L, Cvector, *Cvec ;
76         /* get input L */
77         Lin = cs_dl_mex_get_sparse (&Lmatrix, 1, 1, pargin [0]) ;
78         /* make a copy of L (this can take more work than updating L itself) */
79         lnz = Lin->p [n] ;
80         L = cs_dl_spalloc (n, n, lnz, 1, 0) ;
81         for (j = 0 ; j <= n ; j++) L->p [j] = Lin->p [j] ;
82         for (k = 0 ; k < lnz ; k++) L->i [k] = Lin->i [k] ;
83         for (k = 0 ; k < lnz ; k++) L->x [k] = Lin->x [k] ;
84         cs_mex_check (0, n, -1, 0, 1, 1, pargin [1]) ;     /* get C */
85         C = cs_dl_mex_get_sparse (&Cmatrix, 0, 1, pargin [1]) ;
86 
87         /* do the update one column at a time */
88         Cvec = &Cvector ;
89         Cvec->m = n ;
90         Cvec->n = 1 ;
91         Cvec->p = cp ;
92         Cvec->nz = -1 ;
93         cp [0] = 0 ;
94         for (k = 0 ; k < C->n ; k++)
95         {
96             /* extract C(:,k) */
97             cp [1] = C->p [k+1] - C->p [k] ;
98             Cvec->nzmax = cp [1] ;
99             Cvec->i = C->i + C->p [k] ;
100             Cvec->x = C->x + C->p [k] ;
101             /* update/downdate */
102             ok = cs_dl_updown (L, sigma, Cvec, parent) ;
103             if (!ok) mexErrMsgTxt ("matrix is not positive definite") ;
104         }
105         /* return new L */
106         pargout [0] = cs_dl_mex_put_sparse (&L) ;
107     }
108 }
109