//------------------------------------------------------------------------------ // CCOLAMD/MATLAB/ccolamdtestmex.c: MATLAB test for CCOLAMD //------------------------------------------------------------------------------ // CCOLAMD, Copyright (c) 2005-2022, Univ. of Florida, All Rights Reserved. // Authors: Timothy A. Davis, Sivasankaran Rajamanickam, and Stefan Larimore. // SPDX-License-Identifier: BSD-3-clause //------------------------------------------------------------------------------ /* * This MATLAB mexFunction is for testing only. It is not meant for * production use. See ccolamdmex.c and ccolamd.m instead. * * Usage: * * [ P, stats ] = ccolamdtest (A, knobs) ; * * The stats vector is optional. knobs is required: * * knobs (1) order for LU if nonzero, Cholesky otherwise. default 0 * knobs (2) spumoni, default 0. * knobs (3) dense row control. default 10 * knobs (4) dense column control. default 10 * knobs (5) aggresive absorption if nonzero. default: 1 * * knobs (6) for testing only. Controls the workspace used. * knobs (7) for testing only. Controls how the input matrix is * jumbled prior to calling colamd, to test its error * handling capability. * * see ccolamd.c for a description of the stats array. * */ /* ========================================================================== */ /* === Include files ======================================================== */ /* ========================================================================== */ #include "ccolamd.h" #include "mex.h" #include "matrix.h" #include #include /* Here only for testing */ #undef MIN #undef MAX #define MIN(a,b) (((a) < (b)) ? (a) : (b)) #define MAX(a,b) (((a) > (b)) ? (a) : (b)) #define CCOLAMD_MIN_MEMORY(nnz,n_row,n_col) \ (MAX (2 * nnz, 4 * n_col) + \ 8*n_col + 6*n_row + n_col + (nnz / 5) \ + ((3 * n_col) + 1) + 5 * (n_col + 1) + n_row) /* ========================================================================== */ /* === dump_matrix ========================================================== */ /* ========================================================================== */ static void dump_matrix ( int64_t A [ ], int64_t p [ ], int64_t n_row, int64_t n_col, int64_t Alen, int64_t limit ) { int64_t col, k, row ; mexPrintf ("dump matrix: nrow %d ncol %d Alen %d\n", n_row, n_col, Alen) ; for (col = 0 ; col < MIN (n_col, limit) ; col++) { mexPrintf ("column %d, p[col] %d, p [col+1] %d, length %d\n", col, p [col], p [col+1], p [col+1] - p [col]) ; for (k = p [col] ; k < p [col+1] ; k++) { row = A [k] ; mexPrintf (" %d", row) ; } mexPrintf ("\n") ; } } /* ========================================================================== */ /* === ccolamd mexFunction ================================================== */ /* ========================================================================== */ void mexFunction ( /* === Parameters ======================================================= */ int nargout, /* number of left-hand sides */ mxArray *pargout [ ], /* left-hand side matrices */ int nargin, /* number of right--hand sides */ const mxArray *pargin [ ] /* right-hand side matrices */ ) { /* === Local variables ================================================== */ int64_t *A ; /* ccolamd's copy of the matrix and workspace */ int64_t *p ; /* ccolamd's copy of the column pointers */ int64_t Alen ; /* size of A */ int64_t n_col ; /* number of columns of A */ int64_t n_row ; /* number of rows of A */ int64_t nnz ; /* number of entries in A */ int full ; /* TRUE if input matrix full, FALSE if sparse */ double knobs [CCOLAMD_KNOBS] ; /* ccolamd user-controllable parameters */ double *out_perm ; /* output permutation vector */ double *out_stats ; /* output stats vector */ double *in_knobs ; /* input knobs vector */ int64_t i ; /* loop counter */ mxArray *Ainput ; /* input matrix handle */ int spumoni ; /* verbosity variable */ int64_t stats2 [CCOLAMD_STATS] ; /* stats for ccolamd */ int64_t *cp, *cp_end, col, length ; int64_t *stats ; int result, ok ; stats = stats2 ; /* === Check inputs ===================================================== */ if (nargin != 3 || nargout < 0 || nargout > 2) { mexErrMsgTxt ( "ccolamdtest: incorrect number of input and/or output arguments") ; } /* for testing we require all 7 knobs */ if (mxGetNumberOfElements (pargin [1]) != 7) { mexErrMsgTxt ("ccolamdtest: must have all 7 knobs for testing") ; } /* === Get knobs ======================================================== */ ccolamd_l_set_defaults (knobs) ; spumoni = 0 ; in_knobs = mxGetPr (pargin [1]) ; knobs [CCOLAMD_LU] = (in_knobs [0] != 0) ; knobs [CCOLAMD_DENSE_ROW] = in_knobs [1] ; knobs [CCOLAMD_DENSE_COL] = in_knobs [2] ; knobs [CCOLAMD_AGGRESSIVE] = (in_knobs [3] != 0) ; spumoni = (in_knobs [4] != 0) ; /* print knob settings if spumoni is set */ if (spumoni) { mexPrintf ("\nccolamd version %d.%d, %s:\nknobs(1): %g, order for %s\n", CCOLAMD_MAIN_VERSION, CCOLAMD_SUB_VERSION, CCOLAMD_DATE, in_knobs [0], (knobs [CCOLAMD_LU] != 0) ? "lu(A)" : "chol(A'*A)") ; if (knobs [CCOLAMD_DENSE_ROW] >= 0) { mexPrintf ("knobs(2): %g, rows with > max(16,%g*sqrt(size(A,2)))" " entries removed\n", in_knobs [1], knobs [CCOLAMD_DENSE_ROW]) ; } else { mexPrintf ("knobs(2): %g, no dense rows removed\n", in_knobs [1]) ; } if (knobs [CCOLAMD_DENSE_COL] >= 0) { mexPrintf ("knobs(3): %g, cols with > max(16,%g*sqrt(min(size(A)))" " entries removed\n", in_knobs [2], knobs [CCOLAMD_DENSE_COL]) ; } else { mexPrintf ("knobs(3): no dense columns removed\n", in_knobs [2]) ; } mexPrintf ("knobs(4): %g, aggressive absorption: %s\n", in_knobs [3], (knobs [CCOLAMD_AGGRESSIVE] != 0) ? "yes" : "no") ; mexPrintf ("knobs(5): %g, statistics and knobs printed\n", in_knobs [4]) ; mexPrintf ("Testing: %g %g\n", in_knobs [5], in_knobs[6]) ; } /* === If A is full, convert to a sparse matrix ========================= */ Ainput = (mxArray *) pargin [0] ; if (mxGetNumberOfDimensions (Ainput) != 2) { mexErrMsgTxt ("ccolamd: input matrix must be 2-dimensional") ; } full = !mxIsSparse (Ainput) ; if (full) { mexCallMATLAB (1, &Ainput, 1, (mxArray **) pargin, "sparse") ; } /* === Allocate workspace for ccolamd =================================== */ /* get size of matrix */ n_row = mxGetM (Ainput) ; n_col = mxGetN (Ainput) ; /* get column pointer vector so we can find nnz */ p = (int64_t *) mxCalloc (n_col+1, sizeof (int64_t)) ; (void) memcpy (p, mxGetJc (Ainput), (n_col+1)*sizeof (int64_t)) ; nnz = p [n_col] ; Alen = (int64_t) ccolamd_l_recommended (nnz, n_row, n_col) ; if (Alen == 0) { mexErrMsgTxt ("ccolamd: problem too large") ; } /* === Modify size of Alen if testing =================================== */ /* knobs [5] amount of workspace given to colamd. < 0 : TIGHT memory > 0 : MIN + knob [3] - 1 == 0 : RECOMMENDED memory */ /* get knob [5], if negative */ if (in_knobs [5] < 0) { Alen = CCOLAMD_MIN_MEMORY (nnz, n_row, n_col) + n_col ; } else if (in_knobs [5] > 0) { Alen = CCOLAMD_MIN_MEMORY (nnz, n_row, n_col) + in_knobs [5] - 1 ; } /* otherwise, we use the recommended amount set above */ /* === Copy input matrix into workspace ================================= */ A = (int64_t *) mxCalloc (Alen, sizeof (int64_t)) ; (void) memcpy (A, mxGetIr (Ainput), nnz*sizeof (int64_t)) ; if (full) { mxDestroyArray (Ainput) ; } /* === Jumble matrix ==================================================== */ /* knobs [6] FOR TESTING ONLY: Specifies how to jumble matrix 0 : No jumbling 1 : Make n_row less than zero 2 : Make first pointer non-zero 3 : Make column pointers not non-decreasing 4 : Make a column pointer greater or equal to Alen 5 : Make row indices not strictly increasing 6 : Make a row index greater or equal to n_row 7 : Set A = NULL 8 : Set p = NULL 9 : Repeat row index 10: make row indices not sorted 11: jumble columns massively (note this changes the pattern of the matrix A.) 12: Set stats = NULL 13: Make n_col less than zero */ /* jumble appropriately */ switch ((int) in_knobs [6]) { case 0 : if (spumoni) { mexPrintf ("ccolamdtest: no errors expected\n") ; } result = 1 ; /* no errors */ break ; case 1 : if (spumoni) { mexPrintf ("ccolamdtest: nrow out of range\n") ; } result = 0 ; /* nrow out of range */ n_row = -1 ; break ; case 2 : if (spumoni) { mexPrintf ("ccolamdtest: p [0] nonzero\n") ; } result = 0 ; /* p [0] must be zero */ p [0] = 1 ; break ; case 3 : if (spumoni) { mexPrintf ("colamdtest: negative length last column\n") ; } result = (n_col == 0) ; /* p must be monotonically inc. */ p [n_col] = p [0] ; break ; case 4 : if (spumoni) { mexPrintf ("colamdtest: Alen too small\n") ; } result = 0 ; /* out of memory */ p [n_col] = Alen ; break ; case 5 : if (spumoni) { mexPrintf ("colamdtest: row index out of range (-1)\n") ; } if (nnz > 0) /* row index out of range */ { result = 0 ; A [nnz-1] = -1 ; } else { if (spumoni) { mexPrintf ("Note: no row indices to put out of range\n") ; } result = 1 ; } break ; case 6 : if (spumoni) { mexPrintf ("ccolamdtest: row index out of range (n_row)\n") ; } if (nnz > 0) /* row index out of range */ { if (spumoni) { mexPrintf ("Changing A[nnz-1] from %d to %d\n", A [nnz-1], n_row) ; } result = 0 ; A [nnz-1] = n_row ; } else { if (spumoni) { mexPrintf ("Note: no row indices to put out of range\n") ; } result = 1 ; } break ; case 7 : if (spumoni) { mexPrintf ("ccolamdtest: A not present\n") ; } result = 0 ; /* A not present */ A = (int64_t *) NULL ; break ; case 8 : if (spumoni) { mexPrintf ("ccolamdtest: p not present\n") ; } result = 0 ; /* p not present */ p = (int64_t *) NULL ; break ; case 9 : if (spumoni) { mexPrintf ("ccolamdtest: duplicate row index\n") ; } result = 1 ; /* duplicate row index */ for (col = 0 ; col < n_col ; col++) { length = p [col+1] - p [col] ; if (length > 1) { A [p [col]] = A [p [col] + 1] ; if (spumoni) { mexPrintf ("Made duplicate row %d in col %d\n", A [p [col] + 1], col) ; } break ; } } if (spumoni > 1) { dump_matrix (A, p, n_row, n_col, Alen, col+2) ; } break ; case 10 : if (spumoni) { mexPrintf ("ccolamdtest: unsorted column\n") ; } result = 1 ; /* jumbled columns */ for (col = 0 ; col < n_col ; col++) { length = p [col+1] - p [col] ; if (length > 1) { i = A[p [col]] ; A [p [col]] = A[p [col] + 1] ; A [p [col] + 1] = i ; if (spumoni) { mexPrintf ("Unsorted column %d \n", col) ; } break ; } } if (spumoni > 1) { dump_matrix (A, p, n_row, n_col, Alen, col+2) ; } break ; case 11 : if (spumoni) { mexPrintf ("ccolamdtest: massive jumbling\n") ; } result = 1 ; /* massive jumbling, but no errors */ srand (1) ; for (i = 0 ; i < n_col ; i++) { cp = &A [p [i]] ; cp_end = &A [p [i+1]] ; while (cp < cp_end) { *cp++ = rand() % n_row ; } } if (spumoni > 1) { dump_matrix (A, p, n_row, n_col, Alen, n_col) ; } break ; case 12 : if (spumoni) { mexPrintf ("ccolamdtest: stats not present\n") ; } result = 0 ; /* stats not present */ stats = (int64_t *) NULL ; break ; case 13 : if (spumoni) { mexPrintf ("ccolamdtest: ncol out of range\n") ; } result = 0 ; /* ncol out of range */ n_col = -1 ; break ; } /* === Order the columns (destroys A) =================================== */ ok = ccolamd_l (n_row, n_col, Alen, A, p, knobs, stats, NULL) ; if (spumoni) { ccolamd_l_report (stats) ; } /* === Return the stats vector ========================================== */ if (nargout == 2) { pargout [1] = mxCreateDoubleMatrix (1, CCOLAMD_STATS, mxREAL) ; out_stats = mxGetPr (pargout [1]) ; for (i = 0 ; i < CCOLAMD_STATS ; i++) { out_stats [i] = (stats == NULL) ? (-1) : (stats [i]) ; } /* fix stats (5) and (6), for 1-based information on jumbled matrix. */ /* note that this correction doesn't occur if csymamd returns FALSE */ out_stats [CCOLAMD_INFO1] ++ ; out_stats [CCOLAMD_INFO2] ++ ; } mxFree (A) ; if (ok) { /* === Return the permutation vector ================================ */ pargout [0] = mxCreateDoubleMatrix (1, n_col, mxREAL) ; out_perm = mxGetPr (pargout [0]) ; for (i = 0 ; i < n_col ; i++) { /* ccolamd is 0-based, but MATLAB expects this to be 1-based */ out_perm [i] = p [i] + 1 ; } if (!result) { ccolamd_l_report (stats) ; mexErrMsgTxt ("ccolamd should have returned TRUE\n") ; } } else { /* return p = -1 if ccolamd failed */ pargout [0] = mxCreateDoubleMatrix (1, 1, mxREAL) ; out_perm = mxGetPr (pargout [0]) ; out_perm [0] = -1 ; if (result) { ccolamd_l_report (stats) ; mexErrMsgTxt ("ccolamd should have returned FALSE\n") ; } } mxFree (p) ; }