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