1 /* ========================================================================== */ 2 /* === BTF package ========================================================== */ 3 /* ========================================================================== */ 4 5 /* BTF_MAXTRANS: find a column permutation Q to give A*Q a zero-free diagonal 6 * BTF_STRONGCOMP: find a symmetric permutation P to put P*A*P' into block 7 * upper triangular form. 8 * BTF_ORDER: do both of the above (btf_maxtrans then btf_strongcomp). 9 * 10 * By Tim Davis. Copyright (c) 2004-2007, University of Florida. 11 * with support from Sandia National Laboratories. All Rights Reserved. 12 */ 13 14 15 /* ========================================================================== */ 16 /* === BTF_MAXTRANS ========================================================= */ 17 /* ========================================================================== */ 18 19 /* BTF_MAXTRANS: finds a permutation of the columns of a matrix so that it has a 20 * zero-free diagonal. The input is an m-by-n sparse matrix in compressed 21 * column form. The array Ap of size n+1 gives the starting and ending 22 * positions of the columns in the array Ai. Ap[0] must be zero. The array Ai 23 * contains the row indices of the nonzeros of the matrix A, and is of size 24 * Ap[n]. The row indices of column j are located in Ai[Ap[j] ... Ap[j+1]-1]. 25 * Row indices must be in the range 0 to m-1. Duplicate entries may be present 26 * in any given column. The input matrix is not checked for validity (row 27 * indices out of the range 0 to m-1 will lead to an undeterminate result - 28 * possibly a core dump, for example). Row indices in any given column need 29 * not be in sorted order. However, if they are sorted and the matrix already 30 * has a zero-free diagonal, then the identity permutation is returned. 31 * 32 * The output of btf_maxtrans is an array Match of size n. If row i is matched 33 * with column j, then A(i,j) is nonzero, and then Match[i] = j. If the matrix 34 * is structurally nonsingular, all entries in the Match array are unique, and 35 * Match can be viewed as a column permutation if A is square. That is, column 36 * k of the original matrix becomes column Match[k] of the permuted matrix. In 37 * MATLAB, this can be expressed as (for non-structurally singular matrices): 38 * 39 * Match = maxtrans (A) ; 40 * B = A (:, Match) ; 41 * 42 * except of course here the A matrix and Match vector are all 0-based (rows 43 * and columns in the range 0 to n-1), not 1-based (rows/cols in range 1 to n). 44 * The MATLAB dmperm routine returns a row permutation. See the maxtrans 45 * mexFunction for more details. 46 * 47 * If row i is not matched to any column, then Match[i] is == -1. The 48 * btf_maxtrans routine returns the number of nonzeros on diagonal of the 49 * permuted matrix. 50 * 51 * In the MATLAB mexFunction interface to btf_maxtrans, 1 is added to the Match 52 * array to obtain a 1-based permutation. Thus, in MATLAB where A is m-by-n: 53 * 54 * q = maxtrans (A) ; % has entries in the range 0:n 55 * q % a column permutation (only if sprank(A)==n) 56 * B = A (:, q) ; % permuted matrix (only if sprank(A)==n) 57 * sum (q > 0) ; % same as "sprank (A)" 58 * 59 * This behaviour differs from p = dmperm (A) in MATLAB, which returns the 60 * matching as p(j)=i if row i and column j are matched, and p(j)=0 if column j 61 * is unmatched. 62 * 63 * p = dmperm (A) ; % has entries in the range 0:m 64 * p % a row permutation (only if sprank(A)==m) 65 * B = A (p, :) ; % permuted matrix (only if sprank(A)==m) 66 * sum (p > 0) ; % definition of sprank (A) 67 * 68 * This algorithm is based on the paper "On Algorithms for obtaining a maximum 69 * transversal" by Iain Duff, ACM Trans. Mathematical Software, vol 7, no. 1, 70 * pp. 315-330, and "Algorithm 575: Permutations for a zero-free diagonal", 71 * same issue, pp. 387-390. Algorithm 575 is MC21A in the Harwell Subroutine 72 * Library. This code is not merely a translation of the Fortran code into C. 73 * It is a completely new implementation of the basic underlying method (depth 74 * first search over a subgraph with nodes corresponding to columns matched so 75 * far, and cheap matching). This code was written with minimal observation of 76 * the MC21A/B code itself. See comments below for a comparison between the 77 * maxtrans and MC21A/B codes. 78 * 79 * This routine operates on a column-form matrix and produces a column 80 * permutation. MC21A uses a row-form matrix and produces a row permutation. 81 * The difference is merely one of convention in the comments and interpretation 82 * of the inputs and outputs. If you want a row permutation, simply pass a 83 * compressed-row sparse matrix to this routine and you will get a row 84 * permutation (just like MC21A). Similarly, you can pass a column-oriented 85 * matrix to MC21A and it will happily return a column permutation. 86 */ 87 88 #ifndef _BTF_H 89 #define _BTF_H 90 91 /* make it easy for C++ programs to include BTF */ 92 #ifdef __cplusplus 93 extern "C" { 94 #endif 95 96 #include "SuiteSparse_config.h" 97 98 int btf_maxtrans /* returns # of columns matched */ 99 ( 100 /* --- input, not modified: --- */ 101 int nrow, /* A is nrow-by-ncol in compressed column form */ 102 int ncol, 103 int Ap [ ], /* size ncol+1 */ 104 int Ai [ ], /* size nz = Ap [ncol] */ 105 double maxwork, /* maximum amount of work to do is maxwork*nnz(A); no limit 106 * if <= 0 */ 107 108 /* --- output, not defined on input --- */ 109 double *work, /* work = -1 if maxwork > 0 and the total work performed 110 * reached the maximum of maxwork*nnz(A). 111 * Otherwise, work = the total work performed. */ 112 113 int Match [ ], /* size nrow. Match [i] = j if column j matched to row i 114 * (see above for the singular-matrix case) */ 115 116 /* --- workspace, not defined on input or output --- */ 117 int Work [ ] /* size 5*ncol */ 118 ) ; 119 120 /* long integer version (all "int" parameters become "SuiteSparse_long") */ 121 SuiteSparse_long btf_l_maxtrans (SuiteSparse_long, SuiteSparse_long, 122 SuiteSparse_long *, SuiteSparse_long *, double, double *, 123 SuiteSparse_long *, SuiteSparse_long *) ; 124 125 126 /* ========================================================================== */ 127 /* === BTF_STRONGCOMP ======================================================= */ 128 /* ========================================================================== */ 129 130 /* BTF_STRONGCOMP finds the strongly connected components of a graph, returning 131 * a symmetric permutation. The matrix A must be square, and is provided on 132 * input in compressed-column form (see BTF_MAXTRANS, above). The diagonal of 133 * the input matrix A (or A*Q if Q is provided on input) is ignored. 134 * 135 * If Q is not NULL on input, then the strongly connected components of A*Q are 136 * found. Q may be flagged on input, where Q[k] < 0 denotes a flagged column k. 137 * The permutation is j = BTF_UNFLIP (Q [k]). On output, Q is modified (the 138 * flags are preserved) so that P*A*Q is in block upper triangular form. 139 * 140 * If Q is NULL, then the permutation P is returned so that P*A*P' is in upper 141 * block triangular form. 142 * 143 * The vector R gives the block boundaries, where block b is in rows/columns 144 * R[b] to R[b+1]-1 of the permuted matrix, and where b ranges from 1 to the 145 * number of strongly connected components found. 146 */ 147 148 int btf_strongcomp /* return # of strongly connected components */ 149 ( 150 /* input, not modified: */ 151 int n, /* A is n-by-n in compressed column form */ 152 int Ap [ ], /* size n+1 */ 153 int Ai [ ], /* size nz = Ap [n] */ 154 155 /* optional input, modified (if present) on output: */ 156 int Q [ ], /* size n, input column permutation */ 157 158 /* output, not defined on input */ 159 int P [ ], /* size n. P [k] = j if row and column j are kth row/col 160 * in permuted matrix. */ 161 162 int R [ ], /* size n+1. block b is in rows/cols R[b] ... R[b+1]-1 */ 163 164 /* workspace, not defined on input or output */ 165 int Work [ ] /* size 4n */ 166 ) ; 167 168 SuiteSparse_long btf_l_strongcomp (SuiteSparse_long, SuiteSparse_long *, 169 SuiteSparse_long *, SuiteSparse_long *, SuiteSparse_long *, 170 SuiteSparse_long *, SuiteSparse_long *) ; 171 172 173 /* ========================================================================== */ 174 /* === BTF_ORDER ============================================================ */ 175 /* ========================================================================== */ 176 177 /* BTF_ORDER permutes a square matrix into upper block triangular form. It 178 * does this by first finding a maximum matching (or perhaps a limited matching 179 * if the work is limited), via the btf_maxtrans function. If a complete 180 * matching is not found, BTF_ORDER completes the permutation, but flags the 181 * columns of P*A*Q to denote which columns are not matched. If the matrix is 182 * structurally rank deficient, some of the entries on the diagonal of the 183 * permuted matrix will be zero. BTF_ORDER then calls btf_strongcomp to find 184 * the strongly-connected components. 185 * 186 * On output, P and Q are the row and column permutations, where i = P[k] if 187 * row i of A is the kth row of P*A*Q, and j = BTF_UNFLIP(Q[k]) if column j of 188 * A is the kth column of P*A*Q. If Q[k] < 0, then the (k,k)th entry in P*A*Q 189 * is structurally zero. 190 * 191 * The vector R gives the block boundaries, where block b is in rows/columns 192 * R[b] to R[b+1]-1 of the permuted matrix, and where b ranges from 1 to the 193 * number of strongly connected components found. 194 */ 195 196 int btf_order /* returns number of blocks found */ 197 ( 198 /* --- input, not modified: --- */ 199 int n, /* A is n-by-n in compressed column form */ 200 int Ap [ ], /* size n+1 */ 201 int Ai [ ], /* size nz = Ap [n] */ 202 double maxwork, /* do at most maxwork*nnz(A) work in the maximum 203 * transversal; no limit if <= 0 */ 204 205 /* --- output, not defined on input --- */ 206 double *work, /* return value from btf_maxtrans */ 207 int P [ ], /* size n, row permutation */ 208 int Q [ ], /* size n, column permutation */ 209 int R [ ], /* size n+1. block b is in rows/cols R[b] ... R[b+1]-1 */ 210 int *nmatch, /* # nonzeros on diagonal of P*A*Q */ 211 212 /* --- workspace, not defined on input or output --- */ 213 int Work [ ] /* size 5n */ 214 ) ; 215 216 SuiteSparse_long btf_l_order (SuiteSparse_long, SuiteSparse_long *, 217 SuiteSparse_long *, double , double *, SuiteSparse_long *, 218 SuiteSparse_long *, SuiteSparse_long *, SuiteSparse_long *, 219 SuiteSparse_long *) ; 220 221 222 /* ========================================================================== */ 223 /* === BTF marking of singular columns ====================================== */ 224 /* ========================================================================== */ 225 226 /* BTF_FLIP is a "negation about -1", and is used to mark an integer j 227 * that is normally non-negative. BTF_FLIP (-1) is -1. BTF_FLIP of 228 * a number > -1 is negative, and BTF_FLIP of a number < -1 is positive. 229 * BTF_FLIP (BTF_FLIP (j)) = j for all integers j. UNFLIP (j) acts 230 * like an "absolute value" operation, and is always >= -1. You can test 231 * whether or not an integer j is "flipped" with the BTF_ISFLIPPED (j) 232 * macro. 233 */ 234 235 #define BTF_FLIP(j) (-(j)-2) 236 #define BTF_ISFLIPPED(j) ((j) < -1) 237 #define BTF_UNFLIP(j) ((BTF_ISFLIPPED (j)) ? BTF_FLIP (j) : (j)) 238 239 /* ========================================================================== */ 240 /* === BTF version ========================================================== */ 241 /* ========================================================================== */ 242 243 /* All versions of BTF include these definitions. 244 * As an example, to test if the version you are using is 1.2 or later: 245 * 246 * if (BTF_VERSION >= BTF_VERSION_CODE (1,2)) ... 247 * 248 * This also works during compile-time: 249 * 250 * #if (BTF >= BTF_VERSION_CODE (1,2)) 251 * printf ("This is version 1.2 or later\n") ; 252 * #else 253 * printf ("This is an early version\n") ; 254 * #endif 255 */ 256 257 #define BTF_DATE "May 4, 2016" 258 #define BTF_VERSION_CODE(main,sub) ((main) * 1000 + (sub)) 259 #define BTF_MAIN_VERSION 1 260 #define BTF_SUB_VERSION 2 261 #define BTF_SUBSUB_VERSION 6 262 #define BTF_VERSION BTF_VERSION_CODE(BTF_MAIN_VERSION,BTF_SUB_VERSION) 263 264 #ifdef __cplusplus 265 } 266 #endif 267 #endif 268