|
5164
|
1 /* ========================================================================== */ |
|
|
2 /* === symamd mexFunction =================================================== */ |
|
|
3 /* ========================================================================== */ |
|
|
4 |
|
|
5 /* |
|
|
6 Usage: |
|
|
7 |
|
|
8 P = symamd (A) ; |
|
|
9 |
|
|
10 P = symamd (A, knobs) ; |
|
|
11 |
|
|
12 [ P, stats ] = symamd (A) ; |
|
|
13 |
|
|
14 [ P, stats ] = symamd (A, knobs) ; |
|
|
15 |
|
|
16 Returns a permutation vector P such that the Cholesky factorization of |
|
|
17 A (P,P) tends to be sparser than that of A. This routine provides the same |
|
|
18 functionality as SYMMMD, but tends to be much faster and tends to return a |
|
|
19 better permutation vector. Note that the SYMMMD m-file in |
|
|
20 MATLAB 5.2 also performs a symmetric elimination tree post-ordering. This |
|
|
21 mexFunction does not do this post-ordering. This mexFunction is a |
|
|
22 replacement for the p = sparsfun ('symmmd', A) operation. |
|
|
23 |
|
|
24 A must be square, and is assummed to have a symmetric nonzero pattern. |
|
|
25 Only the nonzero pattern of the lower triangular portion of A is accessed. |
|
|
26 This routine constructs a matrix M such that the nonzero pattern of M'M is |
|
|
27 equal to A (assuming A has symmetric pattern), and then performs a column |
|
|
28 ordering of M using colamd. |
|
|
29 |
|
|
30 The knobs and stats vectors are optional: |
|
|
31 |
|
|
32 knobs (1) rows and columns with more than (knobs(1))*n entries |
|
|
33 are removed prior to ordering, and placed last in |
|
|
34 the output ordering. If knobs is not present, then the |
|
|
35 default of 0.5 is used. |
|
|
36 |
|
|
37 knobs (2) print level, similar to spparms ('spumoni') |
|
|
38 |
|
|
39 stats (1) the number of dense (or empty) rows and columms. These |
|
|
40 are ordered last, in their natural order. |
|
|
41 |
|
|
42 stats (2) (same as stats (1)) |
|
|
43 |
|
|
44 stats (3) the number of garbage collections performed. |
|
|
45 |
|
|
46 stats (4) return status: |
|
|
47 |
|
|
48 0: matrix is a valid MATLAB matrix. |
|
|
49 |
|
|
50 1: matrix has duplicate entries or unsorted columns. |
|
|
51 This should not occur in a valid MATLAB matrix, |
|
|
52 but the ordering proceeded anyway by ignoring the |
|
|
53 duplicate row indices in each column. See |
|
|
54 stats (5:7) for more information. |
|
|
55 |
|
|
56 stats (5) highest numbered column that is unsorted or has |
|
|
57 duplicate entries (zero if none) |
|
|
58 |
|
|
59 stats (6) last seen duplicate or unsorted row index |
|
|
60 (zero if none) |
|
|
61 |
|
|
62 stats (7) number of duplicate or unsorted row indices |
|
|
63 |
|
|
64 Authors: |
|
|
65 |
|
|
66 The authors of the code itself are Stefan I. Larimore and Timothy A. |
|
|
67 Davis (davis@cise.ufl.edu), University of Florida. The algorithm was |
|
|
68 developed in collaboration with John Gilbert, Xerox PARC, and Esmond |
|
|
69 Ng, Oak Ridge National Laboratory. |
|
|
70 |
|
|
71 Date: |
|
|
72 |
|
|
73 September 8, 2003. Version 2.3. |
|
|
74 |
|
|
75 Acknowledgements: |
|
|
76 |
|
|
77 This work was supported by the National Science Foundation, under |
|
|
78 grants DMS-9504974 and DMS-9803599. |
|
|
79 |
|
|
80 Notice: |
|
|
81 |
|
|
82 Copyright (c) 1998-2003 by the University of Florida. |
|
|
83 All Rights Reserved. |
|
|
84 |
|
|
85 See http://www.cise.ufl.edu/research/sparse/colamd (the colamd.c |
|
|
86 file) for the License. |
|
|
87 |
|
|
88 Availability: |
|
|
89 |
|
|
90 The colamd/symamd library is available at |
|
|
91 |
|
|
92 http://www.cise.ufl.edu/research/sparse/colamd/ |
|
|
93 |
|
|
94 This is the http://www.cise.ufl.edu/research/sparse/colamd/symamdmex.c |
|
|
95 file. It requires the colamd.c and colamd.h files. |
|
|
96 |
|
|
97 */ |
|
|
98 |
|
|
99 /* ========================================================================== */ |
|
|
100 /* === Include files ======================================================== */ |
|
|
101 /* ========================================================================== */ |
|
|
102 |
|
|
103 #include "colamd.h" |
|
|
104 #include "mex.h" |
|
|
105 #include "matrix.h" |
|
|
106 #include <stdlib.h> |
|
|
107 |
|
|
108 /* ========================================================================== */ |
|
|
109 /* === symamd mexFunction =================================================== */ |
|
|
110 /* ========================================================================== */ |
|
|
111 |
|
|
112 void mexFunction |
|
|
113 ( |
|
|
114 /* === Parameters ======================================================= */ |
|
|
115 |
|
|
116 int nlhs, /* number of left-hand sides */ |
|
|
117 mxArray *plhs [], /* left-hand side matrices */ |
|
|
118 int nrhs, /* number of right--hand sides */ |
|
|
119 const mxArray *prhs [] /* right-hand side matrices */ |
|
|
120 ) |
|
|
121 { |
|
|
122 /* === Local variables ================================================== */ |
|
|
123 |
|
|
124 int *perm ; /* column ordering of M and ordering of A */ |
|
|
125 int *A ; /* row indices of input matrix A */ |
|
|
126 int *p ; /* column pointers of input matrix A */ |
|
|
127 int n_col ; /* number of columns of A */ |
|
|
128 int n_row ; /* number of rows of A */ |
|
|
129 int full ; /* TRUE if input matrix full, FALSE if sparse */ |
|
|
130 double knobs [COLAMD_KNOBS] ; /* colamd user-controllable parameters */ |
|
|
131 double *out_perm ; /* output permutation vector */ |
|
|
132 double *out_stats ; /* output stats vector */ |
|
|
133 double *in_knobs ; /* input knobs vector */ |
|
|
134 int i ; /* loop counter */ |
|
|
135 mxArray *Ainput ; /* input matrix handle */ |
|
|
136 int spumoni ; /* verbosity variable */ |
|
|
137 int stats [COLAMD_STATS] ; /* stats for symamd */ |
|
|
138 |
|
|
139 /* === Check inputs ===================================================== */ |
|
|
140 |
|
|
141 if (nrhs < 1 || nrhs > 2 || nlhs < 0 || nlhs > 2) |
|
|
142 { |
|
|
143 mexErrMsgTxt ( |
|
|
144 "symamd: incorrect number of input and/or output arguments.") ; |
|
|
145 } |
|
|
146 |
|
|
147 /* === Get knobs ======================================================== */ |
|
|
148 |
|
|
149 colamd_set_defaults (knobs) ; |
|
|
150 spumoni = 0 ; |
|
|
151 |
|
|
152 /* check for user-passed knobs */ |
|
|
153 if (nrhs == 2) |
|
|
154 { |
|
|
155 in_knobs = mxGetPr (prhs [1]) ; |
|
|
156 i = mxGetNumberOfElements (prhs [1]) ; |
|
|
157 if (i > 0) knobs [COLAMD_DENSE_ROW] = in_knobs [COLAMD_DENSE_ROW] ; |
|
|
158 if (i > 1) spumoni = (int) in_knobs [1] ; |
|
|
159 } |
|
|
160 |
|
|
161 /* print knob settings if spumoni > 0 */ |
|
|
162 if (spumoni > 0) |
|
|
163 { |
|
|
164 mexPrintf ("symamd: dense row/col fraction: %f\n", |
|
|
165 knobs [COLAMD_DENSE_ROW]) ; |
|
|
166 } |
|
|
167 |
|
|
168 /* === If A is full, convert to a sparse matrix ========================= */ |
|
|
169 |
|
|
170 Ainput = (mxArray *) prhs [0] ; |
|
|
171 if (mxGetNumberOfDimensions (Ainput) != 2) |
|
|
172 { |
|
|
173 mexErrMsgTxt ("symamd: input matrix must be 2-dimensional.") ; |
|
|
174 } |
|
|
175 full = !mxIsSparse (Ainput) ; |
|
|
176 if (full) |
|
|
177 { |
|
|
178 mexCallMATLAB (1, &Ainput, 1, (mxArray **) prhs, "sparse") ; |
|
|
179 } |
|
|
180 |
|
|
181 /* === Allocate workspace for symamd ==================================== */ |
|
|
182 |
|
|
183 /* get size of matrix */ |
|
|
184 n_row = mxGetM (Ainput) ; |
|
|
185 n_col = mxGetN (Ainput) ; |
|
|
186 if (n_col != n_row) |
|
|
187 { |
|
|
188 mexErrMsgTxt ("symamd: matrix must be square.") ; |
|
|
189 } |
|
|
190 |
|
|
191 A = mxGetIr (Ainput) ; |
|
|
192 p = mxGetJc (Ainput) ; |
|
|
193 perm = (int *) mxCalloc (n_col+1, sizeof (int)) ; |
|
|
194 |
|
|
195 /* === Order the rows and columns of A (does not destroy A) ============= */ |
|
|
196 |
|
|
197 if (!symamd (n_col, A, p, perm, knobs, stats, &mxCalloc, &mxFree)) |
|
|
198 { |
|
|
199 symamd_report (stats) ; |
|
|
200 mexErrMsgTxt ("symamd error!") ; |
|
|
201 } |
|
|
202 |
|
|
203 if (full) |
|
|
204 { |
|
|
205 mxDestroyArray (Ainput) ; |
|
|
206 } |
|
|
207 |
|
|
208 /* === Return the permutation vector ==================================== */ |
|
|
209 |
|
|
210 plhs [0] = mxCreateDoubleMatrix (1, n_col, mxREAL) ; |
|
|
211 out_perm = mxGetPr (plhs [0]) ; |
|
|
212 for (i = 0 ; i < n_col ; i++) |
|
|
213 { |
|
|
214 /* symamd is 0-based, but MATLAB expects this to be 1-based */ |
|
|
215 out_perm [i] = perm [i] + 1 ; |
|
|
216 } |
|
|
217 mxFree (perm) ; |
|
|
218 |
|
|
219 /* === Return the stats vector ========================================== */ |
|
|
220 |
|
|
221 /* print stats if spumoni > 0 */ |
|
|
222 if (spumoni > 0) |
|
|
223 { |
|
|
224 symamd_report (stats) ; |
|
|
225 } |
|
|
226 |
|
|
227 if (nlhs == 2) |
|
|
228 { |
|
|
229 plhs [1] = mxCreateDoubleMatrix (1, COLAMD_STATS, mxREAL) ; |
|
|
230 out_stats = mxGetPr (plhs [1]) ; |
|
|
231 for (i = 0 ; i < COLAMD_STATS ; i++) |
|
|
232 { |
|
|
233 out_stats [i] = stats [i] ; |
|
|
234 } |
|
|
235 |
|
|
236 /* fix stats (5) and (6), for 1-based information on jumbled matrix. */ |
|
|
237 /* note that this correction doesn't occur if symamd returns FALSE */ |
|
|
238 out_stats [COLAMD_INFO1] ++ ; |
|
|
239 out_stats [COLAMD_INFO2] ++ ; |
|
|
240 } |
|
|
241 } |
