csymamdtestmex.c
Go to the documentation of this file.
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 
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 
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 }
CCOLAMD_MAIN_VERSION
#define CCOLAMD_MAIN_VERSION
Definition: ccolamd.h:46
perm
idx_t idx_t idx_t idx_t idx_t * perm
Definition: include/metis.h:223
ccolamd_l_set_defaults
void ccolamd_l_set_defaults(double knobs[CCOLAMD_KNOBS])
col
m col(1)
CCOLAMD_SUB_VERSION
#define CCOLAMD_SUB_VERSION
Definition: ccolamd.h:47
mexFunction
void mexFunction(int nargout, mxArray *pargout[], int nargin, const mxArray *pargin[])
Definition: csymamdtestmex.c:90
result
Values result
Definition: OdometryOptimize.cpp:8
MIN
#define MIN(a, b)
Definition: csymamdtestmex.c:43
A
Matrix< SCALARA, Dynamic, Dynamic, opt_A > A
Definition: bench_gemm.cpp:48
matrix.h
A
Definition: test_numpy_dtypes.cpp:298
CCOLAMD_DENSE_ROW
#define CCOLAMD_DENSE_ROW
Definition: ccolamd.h:63
stats
bool stats
Definition: SolverComparer.cpp:100
Long
#define Long
Definition: csymamdtestmex.c:38
ccolamd.h
CCOLAMD_KNOBS
#define CCOLAMD_KNOBS
Definition: ccolamd.h:57
CCOLAMD_STATS
#define CCOLAMD_STATS
Definition: ccolamd.h:60
CCOLAMD_INFO1
#define CCOLAMD_INFO1
Definition: ccolamd.h:81
return
if n return
Definition: level1_cplx_impl.h:33
row
m row(1)
dump_matrix
static void dump_matrix(Long A[], Long p[], Long n_row, Long n_col, Long Alen, Long limit)
Definition: csymamdtestmex.c:47
CCOLAMD_DATE
#define CCOLAMD_DATE
Definition: ccolamd.h:44
p
float * p
Definition: Tutorial_Map_using.cpp:9
CCOLAMD_INFO2
#define CCOLAMD_INFO2
Definition: ccolamd.h:82
NULL
#define NULL
Definition: ccolamd.c:609
CCOLAMD_DENSE_COL
#define CCOLAMD_DENSE_COL
Definition: ccolamd.h:66
csymamd_l
SuiteSparse_long csymamd_l(SuiteSparse_long n, SuiteSparse_long A[], SuiteSparse_long p[], SuiteSparse_long perm[], double knobs[CCOLAMD_KNOBS], SuiteSparse_long stats[CCOLAMD_STATS], void *(*allocate)(size_t, size_t), void(*release)(void *), SuiteSparse_long cmember[], SuiteSparse_long stype)
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
csymamd_l_report
void csymamd_l_report(SuiteSparse_long stats[CCOLAMD_STATS])
CCOLAMD_AGGRESSIVE
#define CCOLAMD_AGGRESSIVE
Definition: ccolamd.h:69

gtsam
Author(s):
autogenerated on Sat Nov 16 2024 04:02:07