1 /* ========================================================================== */
2 /* === CHOLMOD/MATLAB/resymbol mexFunction ================================== */
3 /* ========================================================================== */
4
5 /* -----------------------------------------------------------------------------
6 * CHOLMOD/MATLAB Module. Copyright (C) 2005-2006, Timothy A. Davis
7 * http://www.suitesparse.com
8 * MATLAB(tm) is a Trademark of The MathWorks, Inc.
9 * -------------------------------------------------------------------------- */
10
11 /* Usage:
12 * L = resymbol (L, A)
13 *
14 * Recompute the symbolic Cholesky factorization of the matrix A. A must be
15 * symmetric. Only tril(A) is used. Entries in L that are not in the Cholesky
16 * factorization of A are removed from L. L can be from an LL' or LDL'
17 * factorization. The numerical values of A are ignored; only its nonzero
18 * pattern is used.
19 */
20
21 /* ========================================================================== */
22
23 #include "cholmod_matlab.h"
24
mexFunction(int nargout,mxArray * pargout[],int nargin,const mxArray * pargin[])25 void mexFunction
26 (
27 int nargout,
28 mxArray *pargout [ ],
29 int nargin,
30 const mxArray *pargin [ ]
31 )
32 {
33 double dummy = 0 ;
34 double *Lx, *Lx2, *Lz, *Lz2 ;
35 Long *Li, *Lp, *Lnz2, *Li2, *Lp2, *ColCount ;
36 cholmod_sparse *A, Amatrix, *Lsparse, *S ;
37 cholmod_factor *L ;
38 cholmod_common Common, *cm ;
39 Long j, s, n, lnz, is_complex ;
40
41 /* ---------------------------------------------------------------------- */
42 /* start CHOLMOD and set parameters */
43 /* ---------------------------------------------------------------------- */
44
45 cm = &Common ;
46 cholmod_l_start (cm) ;
47 sputil_config (SPUMONI, cm) ;
48
49 /* ---------------------------------------------------------------------- */
50 /* check inputs */
51 /* ---------------------------------------------------------------------- */
52
53 if (nargout > 1 || nargin != 2)
54 {
55 mexErrMsgTxt ("usage: L = resymbol (L, A)\n") ;
56 }
57
58 n = mxGetN (pargin [0]) ;
59 if (!mxIsSparse (pargin [0]) || n != mxGetM (pargin [0]))
60 {
61 mexErrMsgTxt ("resymbol: L must be sparse and square") ;
62 }
63 if (n != mxGetM (pargin [1]) || n != mxGetN (pargin [1]))
64 {
65 mexErrMsgTxt ("resymbol: A and L must have same dimensions") ;
66 }
67
68 /* ---------------------------------------------------------------------- */
69 /* get the sparse matrix A */
70 /* ---------------------------------------------------------------------- */
71
72 A = sputil_get_sparse_pattern (pargin [1], &Amatrix, &dummy, cm) ;
73 S = (A == &Amatrix) ? NULL : A ;
74
75 A->stype = -1 ;
76
77 /* A = sputil_get_sparse (pargin [1], &Amatrix, &dummy, -1) ; */
78
79 /* ---------------------------------------------------------------------- */
80 /* construct a copy of the input sparse matrix L */
81 /* ---------------------------------------------------------------------- */
82
83 /* get the MATLAB L */
84 Lp = (Long *) mxGetJc (pargin [0]) ;
85 Li = (Long *) mxGetIr (pargin [0]) ;
86 Lx = mxGetPr (pargin [0]) ;
87 Lz = mxGetPi (pargin [0]) ;
88 is_complex = mxIsComplex (pargin [0]) ;
89
90 /* allocate the CHOLMOD symbolic L */
91 L = cholmod_l_allocate_factor (n, cm) ;
92 L->ordering = CHOLMOD_NATURAL ;
93 ColCount = L->ColCount ;
94 for (j = 0 ; j < n ; j++)
95 {
96 ColCount [j] = Lp [j+1] - Lp [j] ;
97 }
98
99 /* allocate space for a CHOLMOD LDL' packed factor */
100 /* (LL' and LDL' are treated identically) */
101 cholmod_l_change_factor (is_complex ? CHOLMOD_ZOMPLEX : CHOLMOD_REAL,
102 FALSE, FALSE, TRUE, TRUE, L, cm) ;
103
104 /* copy MATLAB L into CHOLMOD L */
105 Lp2 = L->p ;
106 Li2 = L->i ;
107 Lx2 = L->x ;
108 Lz2 = L->z ;
109 Lnz2 = L->nz ;
110 lnz = L->nzmax ;
111 for (j = 0 ; j <= n ; j++)
112 {
113 Lp2 [j] = Lp [j] ;
114 }
115 for (j = 0 ; j < n ; j++)
116 {
117 Lnz2 [j] = Lp [j+1] - Lp [j] ;
118 }
119 for (s = 0 ; s < lnz ; s++)
120 {
121 Li2 [s] = Li [s] ;
122 }
123 for (s = 0 ; s < lnz ; s++)
124 {
125 Lx2 [s] = Lx [s] ;
126 }
127 if (is_complex)
128 {
129 for (s = 0 ; s < lnz ; s++)
130 {
131 Lz2 [s] = Lz [s] ;
132 }
133 }
134
135 /* ---------------------------------------------------------------------- */
136 /* resymbolic factorization */
137 /* ---------------------------------------------------------------------- */
138
139 cholmod_l_resymbol (A, NULL, 0, TRUE, L, cm) ;
140
141 /* ---------------------------------------------------------------------- */
142 /* copy the results back to MATLAB */
143 /* ---------------------------------------------------------------------- */
144
145 Lsparse = cholmod_l_factor_to_sparse (L, cm) ;
146
147 /* return L as a sparse matrix */
148 pargout [0] = sputil_put_sparse (&Lsparse, cm) ;
149
150 /* ---------------------------------------------------------------------- */
151 /* free workspace and the CHOLMOD L, except for what is copied to MATLAB */
152 /* ---------------------------------------------------------------------- */
153
154 cholmod_l_free_factor (&L, cm) ;
155 cholmod_l_free_sparse (&S, cm) ;
156 cholmod_l_finish (cm) ;
157 cholmod_l_print_common (" ", cm) ;
158 /*
159 if (cm->malloc_count != 3 + mxIsComplex (pargout[0])) mexErrMsgTxt ("!") ;
160 */
161 }
162