1 /* ========================================================================== */
2 /* === stongcomp mexFunction ================================================ */
3 /* ========================================================================== */
4
5 /* STRONGCOMP: Find a symmetric permutation to upper block triangular form of
6 * a sparse square matrix.
7 *
8 * Usage:
9 *
10 * [p,r] = strongcomp (A) ;
11 *
12 * [p,q,r] = strongcomp (A,qin) ;
13 *
14 * In the first usage, the permuted matrix is C = A (p,p). In the second usage,
15 * the matrix A (:,qin) is symmetrically permuted to upper block triangular
16 * form, where qin is an input column permutation, and the final permuted
17 * matrix is C = A (p,q). This second usage is equivalent to
18 *
19 * [p,r] = strongcomp (A (:,qin)) ;
20 * q = qin (p) ;
21 *
22 * That is, if qin is not present it is assumed to be qin = 1:n.
23 *
24 * C is the permuted matrix, with a number of blocks equal to length(r)-1.
25 * The kth block is from row/col r(k) to row/col r(k+1)-1 of C.
26 * r(1) is one and the last entry in r is equal to n+1.
27 * The diagonal of A (or A (:,qin)) is ignored.
28 *
29 * strongcomp is normally proceeded by a maximum transversal:
30 *
31 * [p,q,r] = strongcomp (A, maxtrans (A))
32 *
33 * if the matrix has full structural rank. This is identical to
34 *
35 * [p,q,r] = btf (A)
36 *
37 * (except that btf handles the case when A is structurally rank-deficient).
38 * It essentially the same as
39 *
40 * [p,q,r] = dmperm (A)
41 *
42 * except that p, q, and r will differ between btf and dmperm. Both return an
43 * upper block triangular form with a zero-free diagonal. The number and sizes
44 * of the blocks will be identical, but the order of the blocks, and the
45 * ordering within the blocks, can be different. For structurally rank
46 * deficient matrices, dmpmerm returns the maximum matching as a zero-free
47 * diagonal that is above the main diagonal; btf always returns the matching as
48 * the main diagonal (which will thus contain zeros).
49 *
50 * By Tim Davis. Copyright (c) 2004-2007, University of Florida.
51 * with support from Sandia National Laboratories. All Rights Reserved.
52 *
53 * See also maxtrans, btf, dmperm
54 */
55
56 /* ========================================================================== */
57
58 #include "mex.h"
59 #include "btf.h"
60 #define Long SuiteSparse_long
61
mexFunction(int nargout,mxArray * pargout[],int nargin,const mxArray * pargin[])62 void mexFunction
63 (
64 int nargout,
65 mxArray *pargout[],
66 int nargin,
67 const mxArray *pargin[]
68 )
69 {
70 Long b, n, i, k, j, *Ap, *Ai, *P, *R, nblocks, *Work, *Q, jj ;
71 double *Px, *Rx, *Qx ;
72
73 /* ---------------------------------------------------------------------- */
74 /* get inputs and allocate workspace */
75 /* ---------------------------------------------------------------------- */
76
77 if (!((nargin == 1 && nargout <= 2) || (nargin == 2 && nargout <= 3)))
78 {
79 mexErrMsgTxt ("Usage: [p,r] = strongcomp (A)"
80 " or [p,q,r] = strongcomp (A,qin)") ;
81 }
82 n = mxGetM (pargin [0]) ;
83 if (!mxIsSparse (pargin [0]) || n != mxGetN (pargin [0]))
84 {
85 mexErrMsgTxt ("strongcomp: A must be sparse, square, and non-empty") ;
86 }
87
88 /* get sparse matrix A */
89 Ap = (Long *) mxGetJc (pargin [0]) ;
90 Ai = (Long *) mxGetIr (pargin [0]) ;
91
92 /* get output arrays */
93 P = mxMalloc (n * sizeof (Long)) ;
94 R = mxMalloc ((n+1) * sizeof (Long)) ;
95
96 /* get workspace of size 4n (recursive code only needs 2n) */
97 Work = mxMalloc (4*n * sizeof (Long)) ;
98
99 /* get the input column permutation Q */
100 if (nargin == 2)
101 {
102 if (mxGetNumberOfElements (pargin [1]) != n)
103 {
104 mexErrMsgTxt
105 ("strongcomp: qin must be a permutation vector of size n") ;
106 }
107 Qx = mxGetPr (pargin [1]) ;
108 Q = mxMalloc (n * sizeof (Long)) ;
109 /* connvert Qin to 0-based and check validity */
110 for (i = 0 ; i < n ; i++)
111 {
112 Work [i] = 0 ;
113 }
114 for (k = 0 ; k < n ; k++)
115 {
116 j = Qx [k] - 1 ; /* convert to 0-based */
117 jj = BTF_UNFLIP (j) ;
118 if (jj < 0 || jj >= n || Work [jj] == 1)
119 {
120 mexErrMsgTxt
121 ("strongcomp: qin must be a permutation vector of size n") ;
122 }
123 Work [jj] = 1 ;
124 Q [k] = j ;
125 }
126 }
127 else
128 {
129 /* no input column permutation */
130 Q = (Long *) NULL ;
131 }
132
133 /* ---------------------------------------------------------------------- */
134 /* find the strongly-connected components of A */
135 /* ---------------------------------------------------------------------- */
136
137 nblocks = btf_l_strongcomp (n, Ap, Ai, Q, P, R, Work) ;
138
139 /* ---------------------------------------------------------------------- */
140 /* create outputs and free workspace */
141 /* ---------------------------------------------------------------------- */
142
143 /* create P */
144 pargout [0] = mxCreateDoubleMatrix (1, n, mxREAL) ;
145 Px = mxGetPr (pargout [0]) ;
146 for (k = 0 ; k < n ; k++)
147 {
148 Px [k] = P [k] + 1 ; /* convert to 1-based */
149 }
150
151 /* create Q */
152 if (nargin == 2 && nargout > 1)
153 {
154 pargout [1] = mxCreateDoubleMatrix (1, n, mxREAL) ;
155 Qx = mxGetPr (pargout [1]) ;
156 for (k = 0 ; k < n ; k++)
157 {
158 Qx [k] = Q [k] + 1 ; /* convert to 1-based */
159 }
160 }
161
162 /* create R */
163 if (nargout == nargin + 1)
164 {
165 pargout [nargin] = mxCreateDoubleMatrix (1, nblocks+1, mxREAL) ;
166 Rx = mxGetPr (pargout [nargin]) ;
167 for (b = 0 ; b <= nblocks ; b++)
168 {
169 Rx [b] = R [b] + 1 ; /* convert to 1-based */
170 }
171 }
172
173 mxFree (P) ;
174 mxFree (R) ;
175 mxFree (Work) ;
176 if (nargin == 2)
177 {
178 mxFree (Q) ;
179 }
180 }
181