1 /* ========================================================================== */
2 /* === ccolamd mexFunction ================================================== */
3 /* ========================================================================== */
4
5 /* ----------------------------------------------------------------------------
6 * CCOLAMD, Copyright (C), Univ. of Florida. Authors: Timothy A. Davis,
7 * Sivasankaran Rajamanickam, and Stefan Larimore
8 * -------------------------------------------------------------------------- */
9
10 /*
11 * Usage:
12 * p = ccolamd (A) ;
13 * [p stats] = ccolamd (A, knobs, cmember) ;
14 *
15 * See ccolamd.m for a description.
16 */
17
18 /* ========================================================================== */
19 /* === Include files ======================================================== */
20 /* ========================================================================== */
21
22 #include "ccolamd.h"
23 #include "mex.h"
24 #include "matrix.h"
25 #include <stdlib.h>
26 #include <string.h>
27 #define Long SuiteSparse_long
28
29 /* ========================================================================== */
30 /* === ccolamd mexFunction ================================================== */
31 /* ========================================================================== */
32
mexFunction(int nargout,mxArray * pargout[],int nargin,const mxArray * pargin[])33 void mexFunction
34 (
35 /* === Parameters ======================================================= */
36
37 int nargout,
38 mxArray *pargout [ ],
39 int nargin,
40 const mxArray *pargin [ ]
41 )
42 {
43 /* === Local variables ================================================== */
44
45 Long *A ; /* ccolamd's copy of the matrix and workspace */
46 Long *cmember ; /* ccolamd's copy of the constraint set */
47 double *in_cmember ; /* input constraint set */
48 Long *p ; /* ccolamd's copy of the column pointers */
49 Long Alen ; /* size of A */
50 Long cslen ; /* size of CS */
51 Long n_col ; /* number of columns of A */
52 Long n_row ; /* number of rows of A */
53 Long nnz ; /* number of entries in A */
54 Long full ; /* TRUE if input matrix full, FALSE if sparse */
55 double knobs [CCOLAMD_KNOBS] ; /* ccolamd user-controllable parameters */
56 double *out_perm ; /* output permutation vector */
57 double *out_stats ; /* output stats vector */
58 double *in_knobs ; /* input knobs vector */
59 Long i ; /* loop counter */
60 mxArray *Ainput ; /* input matrix handle */
61 Long spumoni ; /* verbosity variable */
62 Long stats [CCOLAMD_STATS] ;/* stats for ccolamd */
63
64 /* === Check inputs ===================================================== */
65
66 if (nargin < 1 || nargin > 3 || nargout < 0 || nargout > 2)
67 {
68 mexErrMsgTxt ("Usage: [p stats] = ccolamd (A, knobs, cmember)") ;
69 }
70
71 /* === Get cmember ====================================================== */
72
73 cmember = NULL ;
74 cslen = 0 ;
75 if (nargin > 2)
76 {
77 in_cmember = mxGetPr (pargin [2]) ;
78 cslen = mxGetNumberOfElements (pargin [2]) ;
79 if (cslen != 0)
80 {
81 cmember = (Long *) mxCalloc (cslen, sizeof (Long)) ;
82 for (i = 0 ; i < cslen ; i++)
83 {
84 /* convert cmember from 1-based to 0-based */
85 cmember[i] = ((Long) in_cmember [i] - 1) ;
86 }
87 }
88 }
89
90 /* === Get knobs ======================================================== */
91
92 ccolamd_l_set_defaults (knobs) ;
93 spumoni = 0 ;
94
95 /* check for user-passed knobs */
96 if (nargin > 1)
97 {
98 in_knobs = mxGetPr (pargin [1]) ;
99 i = mxGetNumberOfElements (pargin [1]) ;
100 if (i > 0) knobs [CCOLAMD_LU] = (in_knobs [0] != 0) ;
101 if (i > 1) knobs [CCOLAMD_DENSE_ROW] = in_knobs [1] ;
102 if (i > 2) knobs [CCOLAMD_DENSE_COL] = in_knobs [2] ;
103 if (i > 3) knobs [CCOLAMD_AGGRESSIVE] = (in_knobs [3] != 0) ;
104 if (i > 4) spumoni = (in_knobs [4] != 0) ;
105 }
106
107 /* print knob settings if spumoni is set */
108 if (spumoni)
109 {
110 mexPrintf ("\nccolamd version %d.%d, %s:\nknobs(1): %g, order for %s\n",
111 CCOLAMD_MAIN_VERSION, CCOLAMD_SUB_VERSION, CCOLAMD_DATE,
112 in_knobs [0],
113 (knobs [CCOLAMD_LU] != 0) ? "lu(A)" : "chol(A'*A)") ;
114 if (knobs [CCOLAMD_DENSE_ROW] >= 0)
115 {
116 mexPrintf ("knobs(2): %g, rows with > max(16,%g*sqrt(size(A,2)))"
117 " entries removed\n", in_knobs [1], knobs [CCOLAMD_DENSE_ROW]) ;
118 }
119 else
120 {
121 mexPrintf ("knobs(2): %g, no dense rows removed\n", in_knobs [1]) ;
122 }
123 if (knobs [CCOLAMD_DENSE_COL] >= 0)
124 {
125 mexPrintf ("knobs(3): %g, cols with > max(16,%g*sqrt(min(size(A)))"
126 " entries removed\n", in_knobs [2], knobs [CCOLAMD_DENSE_COL]) ;
127 }
128 else
129 {
130 mexPrintf ("knobs(3): no dense columns removed\n", in_knobs [2]) ;
131 }
132 mexPrintf ("knobs(4): %g, aggressive absorption: %s\n",
133 in_knobs [3], (knobs [CCOLAMD_AGGRESSIVE] != 0) ? "yes" : "no") ;
134 mexPrintf ("knobs(5): %g, statistics and knobs printed\n",
135 in_knobs [4]) ;
136 }
137
138 /* === If A is full, convert to a sparse matrix ========================= */
139
140 Ainput = (mxArray *) pargin [0] ;
141 if (mxGetNumberOfDimensions (Ainput) != 2)
142 {
143 mexErrMsgTxt ("ccolamd: input matrix must be 2-dimensional") ;
144 }
145 full = !mxIsSparse (Ainput) ;
146 if (full)
147 {
148 mexCallMATLAB (1, &Ainput, 1, (mxArray **) pargin, "sparse") ;
149 }
150
151 /* === Allocate workspace for ccolamd =================================== */
152
153 /* get size of matrix */
154 n_row = mxGetM (Ainput) ;
155 n_col = mxGetN (Ainput) ;
156
157 /* get column pointer vector */
158 p = (Long *) mxCalloc (n_col+1, sizeof (Long)) ;
159 (void) memcpy (p, mxGetJc (Ainput), (n_col+1)*sizeof (Long)) ;
160 nnz = p [n_col] ;
161 Alen = (Long) ccolamd_l_recommended (nnz, n_row, n_col) ;
162 if (Alen == 0)
163 {
164 mexErrMsgTxt ("ccolamd: problem too large") ;
165 }
166
167 /* === Copy input matrix into workspace ================================= */
168
169 A = (Long *) mxCalloc (Alen, sizeof (Long)) ;
170 (void) memcpy (A, mxGetIr (Ainput), nnz*sizeof (Long)) ;
171
172 if (full)
173 {
174 mxDestroyArray (Ainput) ;
175 }
176
177 /* Check constraint set size */
178 if (cmember != NULL && cslen != n_col)
179 {
180 mexErrMsgTxt ("ccolamd: cmember must be of length equal to #cols of A");
181 }
182
183 /* === Order the columns (destroys A) =================================== */
184
185 if (!ccolamd_l (n_row, n_col, Alen, A, p, knobs, stats, cmember))
186 {
187 ccolamd_l_report (stats) ;
188 mexErrMsgTxt ("ccolamd error!") ;
189 }
190 mxFree (A) ;
191 mxFree (cmember) ;
192
193 /* === Return the permutation vector ==================================== */
194
195 pargout [0] = mxCreateDoubleMatrix (1, n_col, mxREAL) ;
196 out_perm = mxGetPr (pargout [0]) ;
197 for (i = 0 ; i < n_col ; i++)
198 {
199 /* ccolamd is 0-based, but MATLAB expects this to be 1-based */
200 out_perm [i] = p [i] + 1 ;
201 }
202 mxFree (p) ;
203
204 /* === Return the stats vector ========================================== */
205
206 /* print stats if spumoni is set */
207 if (spumoni)
208 {
209 ccolamd_l_report (stats) ;
210 }
211
212 if (nargout == 2)
213 {
214 pargout [1] = mxCreateDoubleMatrix (1, CCOLAMD_STATS, mxREAL) ;
215 out_stats = mxGetPr (pargout [1]) ;
216 for (i = 0 ; i < CCOLAMD_STATS ; i++)
217 {
218 out_stats [i] = stats [i] ;
219 }
220
221 /* fix stats (5) and (6), for 1-based information on jumbled matrix. */
222 /* note that this correction doesn't occur if symamd returns FALSE */
223 out_stats [CCOLAMD_INFO1] ++ ;
224 out_stats [CCOLAMD_INFO2] ++ ;
225 }
226 }
227