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

gtsam
Author(s):
autogenerated on Sat Nov 16 2024 04:01:58