1 /* ========================================================================== */
2 /* === csymamdtest 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 * This MATLAB mexFunction is for testing only. It is not meant for
12 * production use. See csymamdmex.c and csymamd.m instead.
13 *
14 * Usage:
15 *
16 * [ P, stats ] = csymamdtest (A, knobs, cmember) ;
17 *
18 * The knobs and stats vectors are optional:
19 *
20 * knobs (1) dense row/col control. default 10
21 * knobs (2) spumoni, default 0.
22 * knobs (3) aggresive absorption if nonzero. default 1
23 *
24 * knobs (4) for testing only. Controls how the input matrix is
25 * jumbled prior to calling colamd, to test its error
26 * handling capability.
27 */
28
29 /* ========================================================================== */
30 /* === Include files ======================================================== */
31 /* ========================================================================== */
32
33 #include "ccolamd.h"
34 #include "mex.h"
35 #include "matrix.h"
36 #include <stdlib.h>
37 #include <string.h>
38 #define Long SuiteSparse_long
39
40 #ifdef MIN
41 #undef MIN
42 #endif
43 #define MIN(a,b) (((a) < (b)) ? (a) : (b))
44
45
dump_matrix(Long A[],Long p[],Long n_row,Long n_col,Long Alen,Long limit)46 static void dump_matrix
47 (
48 Long A [ ],
49 Long p [ ],
50 Long n_row,
51 Long n_col,
52 Long Alen,
53 Long limit
54 )
55 {
56 Long col, k, row ;
57
58 mexPrintf ("dump matrix: nrow %d ncol %d Alen %d\n", n_row, n_col, Alen) ;
59
60 if (!A)
61 {
62 mexPrintf ("A not present\n") ;
63 return ;
64 }
65
66 if (!p)
67 {
68 mexPrintf ("p not present\n") ;
69 return ;
70 }
71
72 for (col = 0 ; col < MIN (n_col, limit) ; col++)
73 {
74 mexPrintf ("column %d, p[col] %d, p [col+1] %d, length %d\n",
75 col, p [col], p [col+1], p [col+1] - p [col]) ;
76 for (k = p [col] ; k < p [col+1] ; k++)
77 {
78 row = A [k] ;
79 mexPrintf (" %d", row) ;
80 }
81 mexPrintf ("\n") ;
82 }
83 }
84
85 /* ========================================================================== */
86 /* === csymamd mexFunction ================================================== */
87 /* ========================================================================== */
88
mexFunction(int nargout,mxArray * pargout[],int nargin,const mxArray * pargin[])89 void mexFunction
90 (
91 /* === Parameters ======================================================= */
92
93 int nargout, /* number of left-hand sides */
94 mxArray *pargout [ ], /* left-hand side matrices */
95 int nargin, /* number of right--hand sides */
96 const mxArray *pargin [ ] /* right-hand side matrices */
97 )
98 {
99 /* === Local variables ================================================== */
100
101 Long *perm ; /* column ordering of M and ordering of A */
102 Long *A ; /* row indices of input matrix A */
103 Long *p ; /* column pointers of input matrix A */
104 Long n_col ; /* number of columns of A */
105 Long n_row ; /* number of rows of A */
106 Long full ; /* TRUE if input matrix full, FALSE if sparse */
107 double knobs [CCOLAMD_KNOBS] ; /* ccolamd user-controllable parameters */
108 double *out_perm ; /* output permutation vector */
109 double *out_stats ; /* output stats vector */
110 double *in_knobs ; /* input knobs vector */
111 Long i ; /* loop counter */
112 mxArray *Ainput ; /* input matrix handle */
113 Long spumoni ; /* verbosity variable */
114 Long stats2 [CCOLAMD_STATS] ;/* stats for csymamd */
115
116 Long *cp, *cp_end, result, nnz, col, length, ok ;
117 Long *stats ;
118 stats = stats2 ;
119
120 /* === Check inputs ===================================================== */
121
122 if (nargin != 3 || nargout > 2)
123 {
124 mexErrMsgTxt (
125 "csymamdtest: incorrect number of input and/or output arguments.") ;
126 }
127 /* for testing we require all 4 knobs */
128 if (mxGetNumberOfElements (pargin [1]) != 4)
129 {
130 mexErrMsgTxt ("csymamdtest: must have 4 knobs for testing") ;
131 }
132
133 /* === Get knobs ======================================================== */
134
135 ccolamd_l_set_defaults (knobs) ;
136 spumoni = 0 ;
137
138 in_knobs = mxGetPr (pargin [1]) ;
139
140 i = mxGetNumberOfElements (pargin [1]) ;
141 knobs [CCOLAMD_DENSE_ROW] = in_knobs [0] ;
142 knobs [CCOLAMD_DENSE_COL] = in_knobs [0] ;
143 knobs [CCOLAMD_AGGRESSIVE] = (in_knobs [1] != 0) ;
144 spumoni = (in_knobs [2] != 0) ;
145
146 /* print knob settings if spumoni is set */
147 if (spumoni)
148 {
149 mexPrintf ("\ncsymamd version %d.%d, %s:\n",
150 CCOLAMD_MAIN_VERSION, CCOLAMD_SUB_VERSION, CCOLAMD_DATE) ;
151 if (knobs [CCOLAMD_DENSE_ROW] >= 0)
152 {
153 mexPrintf ("knobs(1): %g, rows/cols with > "
154 "max(16,%g*sqrt(size(A,2))) entries removed\n",
155 in_knobs [0], knobs [CCOLAMD_DENSE_ROW]) ;
156 }
157 else
158 {
159 mexPrintf ("knobs(1): %g, no dense rows removed\n",
160 in_knobs [0]) ;
161 }
162 mexPrintf ("knobs(2): %g, aggressive absorption: %s\n",
163 in_knobs [1], (knobs [CCOLAMD_AGGRESSIVE] != 0) ? "yes" : "no") ;
164 mexPrintf ("knobs(3): %g, statistics and knobs printed\n",
165 in_knobs [2]) ;
166 mexPrintf ("Testing: %g\n", in_knobs [3]) ;
167 }
168
169 /* === If A is full, convert to a sparse matrix ========================= */
170
171 Ainput = (mxArray *) pargin [0] ;
172 if (mxGetNumberOfDimensions (Ainput) != 2)
173 {
174 mexErrMsgTxt ("csymamd: input matrix must be 2-dimensional.") ;
175 }
176 full = !mxIsSparse (Ainput) ;
177 if (full)
178 {
179 mexCallMATLAB (1, &Ainput, 1, (mxArray **) pargin, "sparse") ;
180 }
181
182 /* === Allocate workspace for csymamd =================================== */
183
184 /* get size of matrix */
185 n_row = mxGetM (Ainput) ;
186 n_col = mxGetN (Ainput) ;
187 if (n_col != n_row)
188 {
189 mexErrMsgTxt ("csymamd: matrix must be square.") ;
190 }
191
192 /* p = mxGetJc (Ainput) ; */
193 p = (Long *) mxCalloc (n_col+1, sizeof (Long)) ;
194 (void) memcpy (p, mxGetJc (Ainput), (n_col+1)*sizeof (Long)) ;
195
196 nnz = p [n_col] ;
197 if (spumoni)
198 {
199 mexPrintf ("csymamdtest: nnz %d\n", nnz) ;
200 }
201
202 /* A = mxGetIr (Ainput) ; */
203 A = (Long *) mxCalloc (nnz+1, sizeof (Long)) ;
204 (void) memcpy (A, mxGetIr (Ainput), nnz*sizeof (Long)) ;
205
206 perm = (Long *) mxCalloc (n_col+1, sizeof (Long)) ;
207
208 /* === Jumble matrix ==================================================== */
209
210
211 /*
212 knobs [4] FOR TESTING ONLY: Specifies how to jumble matrix
213 0 : No jumbling
214 1 : (no errors)
215 2 : Make first pointer non-zero
216 3 : Make column pointers not non-decreasing
217 4 : (no errors)
218 5 : Make row indices not strictly increasing
219 6 : Make a row index greater or equal to n_row
220 7 : Set A = NULL
221 8 : Set p = NULL
222 9 : Repeat row index
223 10: make row indices not sorted
224 11: jumble columns massively (note this changes
225 the pattern of the matrix A.)
226 12: Set stats = NULL
227 13: Make n_col less than zero
228 */
229
230 /* jumble appropriately */
231 switch ((Long) in_knobs [3])
232 {
233
234 case 0 :
235 if (spumoni)
236 {
237 mexPrintf ("csymamdtest: no errors expected\n") ;
238 }
239 result = 1 ; /* no errors */
240 break ;
241
242 case 1 :
243 if (spumoni)
244 {
245 mexPrintf ("csymamdtest: no errors expected (1)\n") ;
246 }
247 result = 1 ;
248 break ;
249
250 case 2 :
251 if (spumoni)
252 {
253 mexPrintf ("csymamdtest: p [0] nonzero\n") ;
254 }
255 result = 0 ; /* p [0] must be zero */
256 p [0] = 1 ;
257 break ;
258
259 case 3 :
260 if (spumoni)
261 {
262 mexPrintf ("csymamdtest: negative length last column\n") ;
263 }
264 result = (n_col == 0) ; /* p must be monotonically inc. */
265 p [n_col] = p [0] ;
266 break ;
267
268 case 4 :
269 if (spumoni)
270 {
271 mexPrintf ("csymamdtest: no errors expected (4)\n") ;
272 }
273 result = 1 ;
274 break ;
275
276 case 5 :
277 if (spumoni)
278 {
279 mexPrintf ("csymamdtest: row index out of range (-1)\n") ;
280 }
281 if (nnz > 0) /* row index out of range */
282 {
283 result = 0 ;
284 A [nnz-1] = -1 ;
285 }
286 else
287 {
288 if (spumoni)
289 {
290 mexPrintf ("Note: no row indices to put out of range\n") ;
291 }
292 result = 1 ;
293 }
294 break ;
295
296 case 6 :
297 if (spumoni)
298 {
299 mexPrintf ("csymamdtest: row index out of range (ncol)\n") ;
300 }
301 if (nnz > 0) /* row index out of range */
302 {
303 result = 0 ;
304 A [nnz-1] = n_col ;
305 }
306 else
307 {
308 if (spumoni)
309 {
310 mexPrintf ("Note: no row indices to put out of range\n") ;
311 }
312 result = 1 ;
313 }
314 break ;
315
316 case 7 :
317 if (spumoni)
318 {
319 mexPrintf ("csymamdtest: A not present\n") ;
320 }
321 result = 0 ; /* A not present */
322 A = (Long *) NULL ;
323 break ;
324
325 case 8 :
326 if (spumoni)
327 {
328 mexPrintf ("csymamdtest: p not present\n") ;
329 }
330 result = 0 ; /* p not present */
331 p = (Long *) NULL ;
332 break ;
333
334 case 9 :
335 if (spumoni)
336 {
337 mexPrintf ("csymamdtest: duplicate row index\n") ;
338 }
339 result = 1 ; /* duplicate row index */
340
341 for (col = 0 ; col < n_col ; col++)
342 {
343 length = p [col+1] - p [col] ;
344 if (length > 1)
345 {
346 A [p [col+1]-2] = A [p [col+1] - 1] ;
347 if (spumoni)
348 {
349 mexPrintf ("Made duplicate row %d in col %d\n",
350 A [p [col+1] - 1], col) ;
351 }
352 break ;
353 }
354 }
355
356 if (spumoni > 1)
357 {
358 dump_matrix (A, p, n_row, n_col, nnz, col+2) ;
359 }
360 break ;
361
362 case 10 :
363 if (spumoni)
364 {
365 mexPrintf ("csymamdtest: unsorted column\n") ;
366 }
367 result = 1 ; /* jumbled columns */
368
369 for (col = 0 ; col < n_col ; col++)
370 {
371 length = p [col+1] - p [col] ;
372 if (length > 1)
373 {
374 i = A[p [col]] ;
375 A [p [col]] = A[p [col] + 1] ;
376 A [p [col] + 1] = i ;
377 if (spumoni)
378 {
379 mexPrintf ("Unsorted column %d \n", col) ;
380 }
381 break ;
382 }
383 }
384
385 if (spumoni > 1)
386 {
387 dump_matrix (A, p, n_row, n_col, nnz, col+2) ;
388 }
389 break ;
390
391 case 11 :
392 if (spumoni)
393 {
394 mexPrintf ("csymamdtest: massive jumbling\n") ;
395 }
396 result = 1 ; /* massive jumbling, but no errors */
397 srand (1) ;
398 for (i = 0 ; i < n_col ; i++)
399 {
400 cp = &A [p [i]] ;
401 cp_end = &A [p [i+1]] ;
402 while (cp < cp_end)
403 {
404 *cp++ = rand() % n_row ;
405 }
406 }
407 if (spumoni > 1)
408 {
409 dump_matrix (A, p, n_row, n_col, nnz, n_col) ;
410 }
411 break ;
412
413 case 12 :
414 if (spumoni)
415 {
416 mexPrintf ("csymamdtest: stats not present\n") ;
417 }
418 result = 0 ; /* stats not present */
419 stats = (Long *) NULL ;
420 break ;
421
422 case 13 :
423 if (spumoni)
424 {
425 mexPrintf ("csymamdtest: ncol out of range\n") ;
426 }
427 result = 0 ; /* ncol out of range */
428 n_col = -1 ;
429 break ;
430
431 }
432
433 /* === Order the rows and columns of A (does not destroy A) ============= */
434
435 ok = csymamd_l (n_col, A, p, perm, knobs, stats, &mxCalloc, &mxFree,
436 NULL, -1) ;
437
438 if (full)
439 {
440 mxDestroyArray (Ainput) ;
441 }
442
443 if (spumoni)
444 {
445 csymamd_l_report (stats) ;
446 }
447
448 /* === Return the stats vector ========================================== */
449
450 if (nargout == 2)
451 {
452 pargout [1] = mxCreateDoubleMatrix (1, CCOLAMD_STATS, mxREAL) ;
453 out_stats = mxGetPr (pargout [1]) ;
454 for (i = 0 ; i < CCOLAMD_STATS ; i++)
455 {
456 out_stats [i] = (stats == NULL) ? (-1) : (stats [i]) ;
457 }
458 /* fix stats (5) and (6), for 1-based information on jumbled matrix. */
459 /* note that this correction doesn't occur if csymamd returns FALSE */
460 out_stats [CCOLAMD_INFO1] ++ ;
461 out_stats [CCOLAMD_INFO2] ++ ;
462 }
463
464 mxFree (A) ;
465
466 if (ok)
467 {
468
469 /* === Return the permutation vector ================================ */
470
471 pargout [0] = mxCreateDoubleMatrix (1, n_col, mxREAL) ;
472 out_perm = mxGetPr (pargout [0]) ;
473 for (i = 0 ; i < n_col ; i++)
474 {
475 /* csymamd is 0-based, but MATLAB expects this to be 1-based */
476 out_perm [i] = perm [i] + 1 ;
477 }
478 if (!result)
479 {
480 csymamd_l_report (stats) ;
481 mexErrMsgTxt ("csymamd should have returned TRUE\n") ;
482 }
483 }
484 else
485 {
486
487 /* return p = -1 if csymamd failed */
488 pargout [0] = mxCreateDoubleMatrix (1, 1, mxREAL) ;
489 out_perm = mxGetPr (pargout [0]) ;
490 out_perm [0] = -1 ;
491 if (result)
492 {
493 csymamd_l_report (stats) ;
494 mexErrMsgTxt ("csymamd should have returned FALSE\n") ;
495 }
496 }
497
498 mxFree (p) ;
499 mxFree (perm) ;
500 }
501